MATLAB求微分方程的解
- 格式：docx
- 大小：239.94 KB
- 文档页数：26
实验二 微分方程求解
一、问题背景与实验目的
实际应用问题通过数学建模所归纳而取得的方程，绝大多数都是微分方程，真正能取得代数方程的机遇很少．另一方面，能够求解的微分方程也是十分有限的，专门是高阶方程和偏微分方程（组）．这就要求咱们必需研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精准解），更要研究微分方程（组）的数值解法（近似解）．
对微分方程（组）的解析解法(精准解)，Matlab 有专门的函数能够用，本实验将作必然的介绍．
本实验将要紧研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法．
二、相关函数（命令）及简介
1．dsolve('equ1','equ2',…)：Matlab 求微分方程的解析解．equ一、equ二、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推．
2．simplify(s)：对表达式 s 利用 maple 的化简规那么进行化简．
例如：
syms x
simplify(sin(x)^2 + cos(x)^2)
ans=1
3．[r,how]=simple(s)：由于 Matlab 提供了多种化简规那么，simple 命令确实是对表达式 s 用各类规那么进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规那么．
例如：
syms x
[r,how]=simple(cos(x)^2-sin(x)^2)
r = cos(2*x)
how = combine 4．[T,Y] = solver(odefun,tspan,y0) 求微分方程的数值解．
说明：
(1) 其中的 solver为命令 ode4五、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一．
(2) odefun 是显式常微分方程：00)(),(ytyytfdtdy
(3) 在积分区间 tspan=],[0ftt上，从0t到ft，用初始条件0y求解．
(4) 要取得问题在其他指按时刻点,210,,ttt上的解，那么令 tspan=
],,,[,210ftttt（要求是单调的）．
(5) 因为没有一种算法能够有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver，关于不同的ODE 问题，采纳不同的Solver．
求解器
Solver ODE类型 特点 说明
ode45 非刚性 单步算法；4、5阶Runge-Kutta方程；累计截断误差达3)(x 大部分场合的首选算法
ode23 非刚性 单步算法；2、3阶Runge-Kutta方程；累计截断误差达3)(x 使用于精度较低的情形
ode113 非刚性 多步法；Adams算法；高低精度均可到6310~10 计算时间比 ode45 短
ode23t 适度刚性 采用梯形算法 适度刚性情形
ode15s 刚性 多步法；Gear's反向数值微分；精度中等 若 ode45 失效时，可尝试使用
ode23s 刚性 单步法；2阶 Rosebrock 算法；低精度 当精度较低时，计算时间比 ode15s 短
ode23tb 刚性 梯形算法；低精度 当精度较低时，计算时间比 ode15s 短
(6) 要专门的是：ode23、ode45 是极为经常使用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的经常使用程序，其中：
ode23 采纳龙格-库塔2 阶算法，用3 阶公式作误差估量来调剂步长，具有低等的精度．
ode45 那么采纳龙格-库塔4 阶算法，用5 阶公式作误差估量来调剂步长，具有中等的精度．
5．ezplot(x,y,[tmin,tmax])：符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为 t 的取值范围．
6．inline()：成立一个内联函数．格式：inline('expr', 'var1',
'var2',…) ，注意括号里的表达式要加引号．
例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)
三、实验内容
1. 几个能够直接用 Matlab 求微分方程精准解的例子：
例1：求解微分方程22xxexydxdy，并加以验证．
求解本问题的Matlab 程序为：
syms x y %line1
y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2
diff(y,x)+2*x*y-x*exp(-x^2) %line3
simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4
说明：
(1) 行line1是用命令概念x,y为符号变量．那个地址能够不写，但为确保正确性，建议写上；
(2) 行line2是用命令求出的微分方程的解：
1/2*exp(-x^2)*x^2+exp(-x^2)*C1
(3) 行line3利用所求得的解．那个地址是将解代入原微分方程，结果应该为0，但那个地址给出： -x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)
(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 说明)(xyy的确是微分方程的解．
例2：求微分方程0'xeyxy在初始条件ey2)1(下的特解，并画出解函数的图形．
求解本问题的 Matlab 程序为：
syms x y
y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')
ezplot(y)
微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab格式)，即xeeyx，解函数的图形如图 1：
-6-4-20246-30-20-1001020304050x1/x exp(x)+1/x exp(1)
图1
例3：求微分方程组035yxdtdyeyxdtdxt在初始条件0|,1|00ttyx下的特解，并画出解函数的图形．
求解本问题的 Matlab 程序为：
syms x y t [x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')
simple(x);
simple(y);
ezplot(x,y,[0,]);axis auto
微分方程的特解(式子专门长)和解函数的图形均略．
2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．
例4：求解微分方程初值问题1)0(2222yxxydxdy的数值解，求解范围为区间[0, ]．
fun=inline('-2*y+2*x^2+2*x','x','y');
[x,y]=ode23(fun,[0,],1);
x';
y';
plot(x,y,'o-')
>> x'
ans =
>> y'
ans =
图形结果为图 2． 00.050.10.150.20.250.30.350.40.450.50.60.650.70.750.80.850.90.951
图2
例 5：求解描述振荡器的经典的 Ver der Pol 微分方程
.7,0)0(',1)0(,0)1(222yyydtdyydtyd
分析：令,,121dtdxxyx则.)1(,1221221xxxdtdxxdtdx
先编写函数文件：
function xprime = verderpol(t,x)
global mu;
xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)];
再编写命令文件：
global mu;
mu = 7;
y0=[1;0]
[t,x] = ode45('verderpol',[0,40],y0);
x1=x(:,1);x2=x(:,2);
plot(t,x1)
图形结果为图3． 0510152025303540-2.5-2-1.5-1-0.500.511.522.5
图3
3. 用 Euler 折线法求解
前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方式．
Euler 折线法求解的大体思想是将微分方程初值问题
00)(),,(yxyyxfdxdy
化成一个代数方程，即差分方程，要紧步骤是用差商hxyhxy)()(替代微商dxdy,于是：
)()),(,()()(00xyyxyxfhxyhxykkkk
记)(,1kkkkxyyhxx，从而)(1hxyykk，那么有
1,,2,1,0).,(,),(1100nkyxhfyyhxxxyykkkkkk
例 6：用 Euler 折线法求解微分方程初值问题
1)0(,22yyxydxdy 的数值解(步长h取，求解范围为区间[0,2]．
解：本问题的差分方程为
1,,2,1,0).2),( ),(,,4.0,1,021100nkyxyyxfyxhfyyhxxhyxkkkkkk（其中：
相应的Matlab 程序见附录 1．
附录 1：
clear
f=sym('y+2*x/y^2');
a=0;
b=2;
h=;
n=(b-a)/h+1;
x=0;
y=1;
szj=[x,y];
for i=1:n-1
y=y+h*subs(f,{'x','y'},{x,y});
x=x+h;
szj=[szj;x,y];
end
szj
plot(szj(:,1),szj(:,2))
