MATLAB求解微分方程

• Matlab命令行输入>>rosslersolve(0.2,2,3) • 求解参数a,b,c的值为0.2,2,3时的rossler方程
参数传递方法二
• • • • • • • • Lorenz=@(t,y,beta,rho,sigma)... [-beta*y(1)+y(2)*y(3) -rho*(y(2)-y(3)) -y(1)*y(2)+sigma*y(2)-y(3)]; ts=[0,100];x0=[0 0 1e-10]; [t,y]=ode45(Lorenz,ts,x0,[],8/3,10,28); plot(y(:,1),y(:,2));figure;plot(y(:,1),y(:,3)); figure;plot(y(:,2),y(:,3));
用Matlab求微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 量值 函数 值 ode45 ode23 ode113 ode15s ode23s 由待解 方程写 成的m文件名 ts=[t0,tf]， t0、tf为自 变量的初 值和终值 函数的 初值
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
d 2x 2 dx 2 1000(1 x ) x 0 例3 dt dt x(0) 2; x' (0) 0
STEP1
令 y1=x，y2=y1’
y1 ' y2 2 (1 y1 ) y2 y1 y2 ' 1000 y (0) 2, y (0) 0 1 2
•
2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5
1.5
1
0.5
0 1 0 -1
0 10 20 30 40 50 60 70 80 90 100
1 -2 -3 -2 0 -1
微分方程求解时参数传递方法一
• 假设上例中参数a,b,c的值在运行时确定，我们可 以这样建立函数文件rosslersolve.m • function rosslersolve(a,b,c) • rossler1=@(t,y,a,b,c)[-y(2)-y(3),y(1)+a*y(2),b+(y(1)c)*y(3)]'; • rossler2=@(t,y)rossler1(t,y,a,b,c); • ts=[0 100];x0=[0 0 0]; • [t,y]=ode45(rossler2,ts,x0); • plot3(y(:,1),y(:,2),y(:,3))
图形结果
2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5
0
500
1000
1500
2000
2500
3ห้องสมุดไป่ตู้00
使用匿名函数，STEP2 和STEP3 可以合并为 • vdp2=@(t,y)[y(2),1000*(1-y(1)^2)*y(2)-y(1)]'; • % vdp2为匿名函数 [T,Y]=ode15s(vdp2,[0,3000],[2 ,0]); plot(T,Y(:,1))
结果：x=tan(t+C1)
d 2x t 2 x x 2 y e dt 2 例2 d y 4 x 3 y 4e t dt
MATLAB命令： [x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)',... 'Dy=4*x+3*y+4*exp(-t)'); 结果：x=-6*t*exp(-t)+C1*exp(-t)+C2*exp((1+6^(1/2))*t)+ C3*exp((1-6^(1/2))*t) y=6*t*exp(-t)-C1*exp(-t)+4*C2*exp((1+6^(1/2))*t)+ 2*C2*exp((1+6^(1/2))*t)*6^(1/2)+4*C3*exp((1-6^(1/2))*t)2*C3*exp((1-6^(1/2))*t)*6^(1/2)+1/2*exp(-t)
求微分方程的解析解
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’) 记号: 在表达微分方程时，用字母D表示求微分，D2、D3 等表示求高阶微分. 任何D后所跟的字母为因变量， 自变量可以指定或由系统规则选定为确省
例1
dx 2 1 x dt
MATLAB命令：dsolve(‘Dx=1+x*x’)
20 15 10 5 0 -5 -10 -15 -20 -10
0
10
20
30
40
50
30
20
10
0
-10
-20
-30 -10
0
10
20
30
40
50
30
20
10
0
-10
-20
-30 -20
-15
-10
-5
0
5
10
15
20
化为一阶微分方程组：
• STEP2 建立M文件 • function dy=vdp(t,y) dy=[y(2); 1000*(1-y(1)^2)*y(2)-y(1);]
STEP3 调用MATLAB 函数ODE15S [T,Y]=ode15s('vdp',[0 3000],[2 0]); plot(T,Y(:,1))
dx y z dt dy 例（ 4 Rossler方程） x ay dt dz dt b ( x c) z
选定a=0.3，b=2,c=3 初值x(0)=0,y(0)=0,z(0)=0
• a=0.3;b=2;c=3; • rossler=@(t,y)[-y(2)-y(3),y(1)+a*y(2),b+(y(1)c)*y(3)]'; • ts=[0 100];x0=[0 0 0]; • [t,y]=ode45(rossler,ts,x0); • plot(t,y(:,1),'r',t,y(:,2),'b',t,y(:,3),'g'); • figure • plot3(y(:,1),y(:,2),y(:,3))