MATLAB求解微分方程(实验1)

1、建立M文件函数 function xdot = fun(t,x,y) xdot = [x2(t)；x3(t)；…；f(t, x1(t), x2(t),…xn(t))]； 2、数值计算（执行以下命令）
[t,x1,x2,…,xn]=ode45(‘fun',[t0,tf], [x1(0),x2(0),…,xn(0)])
蓝色曲线 ——y（1）； （原方程解） 红色曲线 ——y（2）；
Van Der Pol Solution 3 2
y(1)and y(2)
选择一组状态变量
x1 y, x2 y,, xn y(n1)
x1 x2 , x2 x3 , xn f (t , x1 , x2 ,, xn )
15
注意
x1 x2 , x2 x3 , xn f (t , x1 , x2 ,, xn )
16
例7
d 2 x dx 2 1000 1 x 2 ) ( x0 dt dt x (0) 2; x ' (0) 0
则微分方程变为一阶微分方程组：
解: 令 y1= x，y2= y1’= x’
y1 ' y2 2 ( y2 ' 1000 1 y1 ) y2 y1 y (0) 2, y (0) 0 2 1
xi 1 xi
f (t , y (t ))dt
xi 1 xi [ f ( xi , y ( xi )) f ( xi 1 , y ( xi 1 ))] 2
故有公式：
h y i 1 y i [ f ( xi , y i ) f ( xi 1 , y i 1 )] 2 y 0 y ( x0 )
ts=[t0,tf]， t0、tf为自 变量的初 值和终值
函数的 初值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法 用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
输入：y=dsolve ('Dy=1+y^2') y1=dsolve('Dy=1+y^2','y(0)=1','x') 输出：y= tan(t-C1) （通解） y1= tan(x+1/4*pi) （特解）
例2
常系数的二阶微分方程
y' '2 y'3 y 0, y(0) 1, y' (0) 0
记号: 在表达微分方程时，用字母 D 表示求微分，D2、D3 等 表示求高阶微分。任何 D 后所跟的字母为因变量，自变量可以 指定或由系统规则选定为缺省。
注意：① y '
‘t’。
Dy，y ' '
D2y
② 自变量名可以省略，默认变量名
3
MATLAB软件求解
例1
dy 1 y2, dx
y (wk.baidu.com0) 1
输入： y=dsolve('D2y-2*Dy-3*y=0','x')
y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')
结果：
y = C1*exp(-x)+C2*exp(3*x)
y = 3/4*exp(-x)+1/4*exp(3*x)
5
例3
非常系数的二阶微分方程
x ( 0) 0 y ( 0) 1
[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y') [x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1') 输出： x =-exp(3*t)*(C1*cos(4*t)-C2*sin(4*t)) y =exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) x =exp(3*t)*sin(4*t) y =exp(3*t)*cos(4*t)
y' f ( x, y ) 对常微分方程： ， 其数值解是指由初始点 x0 开始 y ( x 0 ) y0 的若干离散的 x值处，即对 x0 x1 x 2 x n， 求出准确值 y ( x1 ), y ( x 2 ),, y ( x n ) 的相应近似值 y1 , y2 , , yn。
[t,y]=ode23('vdpol',[0,20],[3,0]); y1=y(:,1); % 原方程的解 y2=y(:,2); plot(t,y1,t,y2, '--' ) % y1(t),y2(t) 曲线图 pause, plot(y1,y2),grid % 相轨迹图，即y2(y1)曲线
20
计算结果
然后继续下一步 i 2 的计算。 y
此即改进的欧拉法。
12
3、使用泰勒公式 以此方法为基础，有龙格-库塔法、线性多步法等方 法。 4、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时 （k为正整数，h为步长），称它是一个k阶公式。 k越大，则数值公式的精度越高。
•欧拉法是一阶公式，改进的欧拉法是二阶公式。
2 1.5 1 0.5 0 -0.5 -1
1、建立m-文件vdp1000.m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
-1.5 -2 -2.5 0
500
1000
1500
实际应用时，与欧拉公式结合使用：
0 y i(1) y i hf ( xi , y i ) h ( k 1) k y i 1 y i [ f ( xi , y i ) f ( xi 1 , y i(1) )] k 0,1,2, 2
k k （ 对于已给的精确度 ， 当满足 yi(11) yi(1) 时， y i 1 yi k11） 取 ，
y 3 ( 0) 1
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12
解 1、建立m-文件rigid.m如下：
function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
x' ' (t ) (1 x(t ) ) x' (t ) x(t ) 0 x(0) 3, x' (0) 0
2
该方程无解析解！
令 y1=x (t), y2 = x’(t)
y1 ' y2 2 y2 ' (1 y1 ) y 2 y1 y1 (0) 3 y 2 (0) 0
8
例 6 求微分方程组的通解.
dx dt 2 x 3 y 3 z dy 4 x 5 y 3z dt dz 4 x 4 y 2 z dt
解 输入命令 ： [x,y,z]=dsolve('Dx=2*x-3*y+3*z', ... 'Dy=4*x-5*y+3*z',... 'Dz=4*x-4*y+2*z', 't'); x=simple(x) % 将x简化 y=simple(y) z=simple(z)
•龙格-库塔法有二阶公式和四阶公式。 •线性多步法有四阶阿达姆斯外插公式和内插公式。 返 回
13
（三）用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 量值 函数 值
ode45 ode23 ode113 ode15s ode23s
由待解 方程写 成的m文件名
x(0) 3, x' (0) 0
x' ' (t ) (1 x 2 (t ))x' (t ) x(t ) 0,
x=dsolve('D2x-(1-x^2)*Dx+x=0', ' x(0)=3,Dx(0)=0')
无解析表达式！
6
MATLAB软件求解 例4 非线性微分方程
x' (t ) x(t ) 1, x(0) 0
2000
2500
3000
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-') 3、结果如图
17
例 8 解微分方程组.
y1 ' y2 y3 y ' y y 2 1 3 y3 ' 0.51y1 y 2 y1 (0) 0, y2 (0) 1,
故有公式：
yi 1 yi hf ( xi , yi ) i 0,1,2, , n - 1 y 0 y ( x0 )
此即欧拉法。
11
2、使用数值积分 对方程y’=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有：
y ( xi 1 ) y ( xi )
实验
Experiments in Mathematics
微 分 方 程 求 解
1
实验目的
1、学会用Matlab求简单微分方程的解析解. 2、学会用Matlab求微分方程的数值解.
实验内容
1、求简单微分方程的解析解. 2、求微分方程的数值解.
实验软件
MATLAB
2
微分方程的解析解
求微分方程（组）的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') 3、结果如图 图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线. 18
例9
Van der pol 方程：
19
（1）编写M文件 ( 文件名为 vdpol.m):
function dy = vdpol(t,y); dy=zeros(2,1); dy(1)=y(2); dy(2)=(1-y(1)^2)*y(2)-y(1); % 或 dy=[y(2);(1-y(1)^2)*y(2)-y(1)];
（2）编写程序如下:（vdj.m）
返 回
10
（二）建立数值解法的一些途径
设 x i 1 x i h, i 0,1,2, n 1, 可用以下离散化方法求 解微分方程： y' f ( x,y) y ( x 0 ) y0
1、用差商代替导数 若步长h较小，则有
y ( x h) y ( x ) y ' ( x) h
结 果 为：x =C3*exp(2*t)+exp(-t)*C1 y =C2*exp(-2*t)+C3*exp(2*t)+exp(-t)*C1 z =C2*exp(-2*t)+C3*exp(2*t)
9
微分方程的数值解
（一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多 得不出一般解。而在实际上对初值问题，一般是要求得 到解在若干个点上满足规定精确度的近似值，或者得到 一个满足精确度要求的便于计算的表达式。 因此，研究常微分方程的数值解法是十分必要的。
14
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量， m-文件中的待解方程组应以x的分量形式写成。 2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组。
y ( n ) f (t , y, y,, y ( n 1) ) y (0), y (0),, y ( n 1) (0)
2 2
x=dsolve('(Dx)^2+x^2=1','x(0)=0')
x = sin(t) -sin(t) 若欲求解的某个数值解，如何求解？
t=pi/2;
eval(x)
7
MATLAB软件求解 例5
输入：
dx dt 3x 4 y dy 4 x 3 y dt