matlab subs函数_[基础向]MATLAB常微分方程的解法

点击蓝字 关注我们

MATLAB

常微分方程的解法

微分方程是模拟自然、生物、工程等现象的重要数学工具。国赛中也有很多与微分方程有关的题目，如2018A高温作业专用服装设计、2014A嫦娥三号软着陆轨道设计与控制策略、2010A储油罐的变位识别与罐容表标定等，都是求解微分方程的定解问题。

★ 符号解法 ★

一、命令详解

 ● y=dsolve(‘equation’, ’cond1, cond2, ...‘, ‘var’);

   求解常微分方程equation满足初值条件cond1, cond2, … 的解。

 ● S=dsolve(‘equation1’, ‘equation2’, …, ‘cond1’, ‘cond2’, …);

  求解多个常微分方程equation1, eqution2, … 满足初值条件cond1, cond2, … 的解，并以结构体的形式输出。

注：常微分方程equation中，符号D表示对变量进行求导，D2，D3，…分别为二阶、三阶导数。

二、案例分析

例1：求解两点边值问题：

代码：

>> y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x');

y =

(31*x^4)/468 - x^3/3 + 125/468

★ 数值解法简介 ★

对于一些微分方程的定解问题，其解析解是不存在的或是难以求得的，因此，我们通常采用一些方法求得微分方程的数值解。所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值。

微分方程和定解条件一起组成定解问题。定解条件通常有两种给法：

(1)初值问题(IVP)：给出积分曲线初始时刻的状态；

(2)边值问题(BVP)：给出积分曲线首末两端的状态。

对于一阶方程，往往只需要初值条件就可以得到方程的解；对于二阶或者二阶以上的微分方程，则可能需要边值条件。

★ 初值问题 ★

一、简介

一阶常微分方程的初值问题：

MATLAB提供了多个一阶常微分方程(组)初值问题数值解的函数，一般函数调用格式为：

 ● [t,y]=solver(fname,tspan,y0, [options]);

其中，t 和 y 分别给出了时间向量和状态向量；solver 为求解常微分方程数值解的函数，细节见下表；fname 为微分方程函数名；tspan 为指定的积分区间；y0 用于指定初值；option 用于改变计算中积分的特性(本文不详述)

二、案例分析

例2：一阶常微分初值问题：

 ● 编写微分方程函数文件 fun1.m

function fy = fun1(x, y)

fy = y.*tan(x)+sec(x);

end

 ● 求解微分方程

[x,y]=ode45(@fun1,[0, 1],pi/2);%求数值解

fy=dsolve('Dy=y*tan(x)+sec(x)','y(0)=pi/2','x');%求精确解

yy=double(subs(fy,'x',x));

plot(x, y,'ro',x,yy,'b-');

 ● 结果

例3：二阶常微分方程初值问题：

● 由于 ode23 和 ode45 都是对一阶常微分方程设计的，对于高阶微分方程，需要先将其转化为一阶常微分方程组，即状态转移方程。令 x1=y,x2=y'，则原方程的转化为：

● 编写微分方程函数文件fun2.m

function fy = fun2(x, t)

fy(1)=x(2);

fy(2)=-x(1)+1-t^2/(2*pi);

end

● 求解微分方程组

[x,y]=ode45(@fun2,[0,3*pi],[0,0]);%求数值解

fy=dsolve('D2y+y=1-t^2/(2*pi)','y(0)=0,Dy(0)=0','t');%求精确解

yy=double(subs(fy,'t',x));

plot(x,y(:,1),'ro',x,yy,'b-');%第一列为所求数值解

● 结果

★ 边值问题 ★

一、命令详解

一阶常微分方程的初值问题：

需要两个特定的条件，除了初值条件外，也可以通过给定边值条件来确定。边值条件一般有三种给定方法(这里不详细给出)。

二、基本命令

 ● solinit = bvpinit(x, yinit)

生成 BVP 求解器的所必须的初始估计值返回解结构体的初始估计值。x 指定边界区间[a,b]上的初始网格，要求x(1)=a，x(end)=b，点格要求单调；yinit是对解的初始猜测，指定为向量或函数。

 ● sol = bvp4c(odefun,bcfun,solinit)

求常微分方程的边界值问题的近似解，返回解结构体。odefun为要求解的函数，指定为函数句柄。bcfun为边界条件，指定为计算边界条件中残差的函数句柄。solinit为解的初始估计值，由函数bvpinit创建。

 ● y = deval(x, sol)

计算 x 中包含的点处的微分方程问题的解 sol。

三、案例分析

例4：考虑二阶常微分方程边值问题：

 ● 将高阶微分方程转化为一阶常微分方程组。令 y1=y,y2=y'，则原方程的转化为：

 ● 编写微分方程函数文件 odefun.m

function fy=odefun(x,y)

fy(1)=y(2);

fy(2)=y(1)+10*x;

end

 ● 编写边界条件函数文件 bcfun.m

function bc=bcfun(ya,yb)

bc(1)=ya(1);

bc(2)=yb(1);

end

 ● 求解微分方程组

solinit = bvpinit(0:0.1:1, [0, 1]);

sol=bvp4c(@odefun, @bcfun, solinit);

x=0:0.02:1;

y=deval(sol,x);  %求数值解

fy=dsolve('D2y-y=10*x','y(0)=0,y(1)=0','x');%求解析解

yy = double(subs(fy, 'x', x));

plot(x, y(1, :), 'ro', x, yy, 'b-');

 ● 结果

- CUMTBSXJM -

本期编辑丨龙临基

责任编辑丨龙临基

责任主编丨张   磊

The End