常微分方程的解法 (四): Matlab 解法
7.1.1 非刚性常微分方程的解法Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如 ode45，ode23， ode113，其中 ode45 采用四五阶 RK 方法，是解非刚性常微分方程的首选方法，ode23 采用二三阶 RK 方法，ode113 采用的是多步法，效率一般比 ode45 高。 Matlab 的工具箱中没有 Euler 方法的功能函数。
（I）对简单的一阶方程的初值问题
https://img-blog.csdnimg.cn/20190430100147247.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
我们自己编写改进的 Euler 方法函数 eulerpro.m 如下：
function =eulerpro(fun,x0,xfinal,y0,n);
if nargin<5,n=50;end
h=(xfinal-x0)/n;
x(1)=x0;y(1)=y0;
for i=1:n
x(i+1)=x(i)+h;
y1=y(i)+h*feval(fun,x(i),y(i));
y2=y(i)+h*feval(fun,x(i+1),y1);
y(i+1)=(y1+y2)/2;
end
例1 用改进的Euler方法求解
https://img-blog.csdnimg.cn/20190430100439679.png
解 编写函数文件 doty.m 如下：
function f=doty(x,y);
f=-2*y+2*x^2+2*x;
在Matlab命令窗口输入：
=eulerpro('doty',0,0.5,1,10)
即可求得数值解。（II）ode23，ode45，ode113的使用 Matlab的函数形式是
=solver('F',tspan，y0)
这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程 y'= f (x, y) 右端的函数。
tspan=
是求解区间，y0是初值。例2 用RK方法求解https://img-blog.csdnimg.cn/20190430100551272.png解 同样地编写函数文件 doty.m 如下：
function f=doty(x,y);
f=-2*y+2*x^2+2*x;
在Matlab命令窗口输入：
=ode45('doty',0,0.5,1)
即可求得数值解。
7.1.2 刚性常微分方程的解法
Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s， ode23t，ode23tb，这些函数的使用同上述非刚性微分方程的功能函数。
7.1.3 高阶微分方程的解法
https://img-blog.csdnimg.cn/20190430101401655.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
（ii）把一阶方程组写成接受两个参数t 和 y ，返回一个列向量的 M 文件 F.m：
function dy=F(t,y);
dy=;
注意：尽管不一定用到参数t 和 y ，M—文件必须接受此两参数。这里向量 dy 必须是列 向量。（iii）用 Matlab 解决此问题的函数形式为
=solver('F',tspan,y0)
这里 solver 为 ode45、ode23、ode113，输入参数 F 是用 M 文件定义的常微分方程组， tspan=是求解区间，y0 是初值列向量。在 Matlab 命令窗口输入 =ode45('F',,) 就得到上述常微分方程的数值解。这里 Y 和时刻 T 是一 一对应的，Y(:,1)是初值问题的 解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。
例 4 求 van der Pol 方程
https://img-blog.csdnimg.cn/20190430101623569.png
的数值解，这里 μ > 0是一参数。
https://img-blog.csdnimg.cn/20190430101709568.png
（ii）书写 M 文件（对于 μ =1）vdp1.m:
function dy=vdp1(t,y);
dy=;
（iii）调用 Matlab 函数。对于初值 y(0) = 2, y'(0) = 0 ，解为
=ode45('vdp1',,);
（iv）观察结果。利用图形输出解的结果；
plot(T,Y(:,1),'-',T,Y(:,2),'--')
title('Solution of van der Pol Equation,mu=1');
xlabel('time t');
ylabel('solution y');
legend('y1','y2');
https://img-blog.csdnimg.cn/20190430101906488.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
例 5 van der Pol 方程， μ =1000 （刚性）解 (i)书写 M 文件 vdp1000.m:
function dy=vdp1000(t,y);
dy=;
（ii）观察结果
=ode15s('vdp1000',,);
plot(t,y(:,1),'o')
title('Solution of van der Pol Equation,mu=1000');
xlabel('time t');
ylabel('solution y(:,1)');
7.2 常微分方程的解析解
在 Matlab 中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令 dsolve。常微分方程在 Matlab 中按如下规定重新表达： 符号 D 表示对变量的求导。Dy 表示对变量 y 求一阶导数，当需要求变量的 n 阶导 数时，用 Dn 表示，D4y 表示对变量 y 求 4 阶导数。 由此，常微分方程 y' '+2y'= y 在 Matlab 中，将写成
D2y+2*Dy=y
7.2.1 求解常微分方程的通解无初边值条件的常微分方程的解就是该方程的通解。其使用格式为：dsolve('diff_equation') dsolve(' diff_equation'，'var')
式中 diff_equation 为待解的常微分方程，第 1 种格式将以变量 t 为自变量进行求解， 第 2 种格式则需定义自变量 var。例 6 试解常微分方程https://img-blog.csdnimg.cn/20190430102152611.png解 编写程序如下：
syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0';
dsolve(diff_equ,'x')
7.2.2 求解常微分方程的初边值问题求解带有初边值条件的常微分方程的使用格式为：dsolve('diff_equation'，'condition1，condition2，…'，'var')
其中 condition1，condition2，… 即为微分方程的初边值条件。例 7 试求微分方程https://img-blog.csdnimg.cn/20190430102316376.png 的解。解 编写程序如下：
y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x')
7.2.3 求解常微分方程组求解常微分方程组的命令格式为：dsolve('diff_equ1，diff_equ2，…'，'var')
dsolve('diff_equ1，diff_equ2，…'，'condition1，condition2，…'，'var')
第 1 种格式用于求解方程组的通解，第 2 种格式可以加上初边值条件，用于具体求解。例 8 试求常微分方程组：https://img-blog.csdnimg.cn/20190430102434331.png的通解和在初边值条件为 f '(2) = 0, f (3) = 3, g(5) = 1的解。解 编写程序如下：
clc,clear
equ1='D2f+3*g=sin(x)';
equ2='Dg+Df=cos(x)';
=dsolve(equ1,equ2,'x')
=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')
7.2.4 求解线性常微分方程组（i）一阶齐次线性微分方程组https://img-blog.csdnimg.cn/2019043010254622.png例 9 试解初值问题https://img-blog.csdnimg.cn/20190430102616792.png解 编写程序如下：
syms t
a=;
x0=;
x=expm(a*t)*x0
（ii）非齐次线性方程组https://img-blog.csdnimg.cn/20190430102716225.png例 10 试解初值问题https://img-blog.csdnimg.cn/20190430102849860.png解 编写程序如下：
clc,clear
syms t s
a=;ft=;
x0=;
x=expm(a*t)*x0+int(expm(a*(t-s))*subs(ft,s),s,0,t);
x=simple(x)
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89703911
页:
[1]