matlab常微分方程的数值解.pptVIP

一、解法步骤： ? 将高阶方程转化为一阶方程组。 （两种情况） ? 建立相应的函数文件。 ? 调用求解函数。 [t,z]=odeij(‘dzdtk,H,z0,tol) 刚性：设有一阶常系数线性微分方程组 y=Ay+f ， 如果它的 Jacobian 矩阵的特征值相差十分悬 殊。 odeij 问题类型 精度 适用对象 ode45 非刚性 中等 多数情况下可优先选用，但不能用来 解刚性问题 ode23 非刚性 较低 可用来解中等刚性问题，或误差允许 范围较宽的问题 ode113 非刚性 低到高 不能解刚性问题，当误差容限要求严 格时效果较 ode45 好 ode15s 刚性 低到高 可用于解刚性问题，当采用 ode45 失败 或效果很差时，可考虑使用 ode23s 刚性 低 可用于解刚性问题，当误差容限较宽 时效果比 ode15s 好 ode23t 适度刚性 低 可用于解刚性问题，但要求无数值衰 减 ode23tb 刚性 低 可用于解刚性问题，当误差容限较宽 时效果比 ode15s 好 二、例题 例 1 在区间 H=[0.1,30] 上满足初值条件 x=0.1 时， 的特解 （ 1 ）转化原方程为一阶方程组。 （ 2 ）建立描述上面微分方程组的函数文件，并设文件名为 dzdt01.m 。 function dz=f(t,z) dz(1)=z(2); dz(2)=z(3); dz(3)=(1/t)*z(3)-(2/t^2)*z(2)+(2/t^3)*z(1)+10*t^3*sin(t); dz=[dz(1);dz(2);dz(3)] x x y y x y x y x sin 10 2 2 6 2 3 ? ? ? ? ? ? ? ? ? ? 0 , 0 , 0 ? ? ? ? ? ? y y y （ 3 ）调用求解。 H=[0.1,30];z0=[0;0;0]; [t,z]=ode45(dzdt01,H,z0); plot(t,z(:,1),m-) plot(t,z(:,1),m-) 0 5 10 15 20 25 30 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 x 10 5 例 2 H=[0,50] （ 1 ）转化方程。 （ 2 ）建立文件，设文件名为 dzdt02.m 。 function dz=f(t,z) c=1; dz(1)=z(2); dz(2)=c*(1-z(1)^2)*z(2)-z(1); dz=[dz(1);dz(2)]; （ 3 ）调用求解。 H=[0,50];z0=[2;1]; [t,z]=ode45(dzdt02,H,z0); plot(t,z(:,1),m-,t,z(:,2),g-.) -3 -2 -1 0 1 2 3 0 ) 1 ( 2 ? ? ? ? ? ? ? y y y c y 1 , 2 , 0 ? ? ? ? y y t 例 3 H=[0,100] （ 1 ）转化方程。 （ 2 ）建立文件，设文件名为 dzdt03.m 。 function dz=f(t,z) a=-0.2; dz(1)=z(2); dz(2)=a*z(2)-sin(t); dz=[dz(1);dz(2)]; （ 3 ）调用求解。 H=[0,100];z0=[0;2]; [t,z]=ode45(dzdt03,H,z0); plot(z(:,1),z(:,2),m-) 0 1 2 3 4 5 6 -1 -0.5 0 0.5 1 1.5 2 t y y y x sin 2 . 0 , ? ? ? ? ? ? 2 , 0 , 0 ? ? ? y x t 例 4 H=[-30 ， 30] （ 1 ）转化方程。 （ 2 ）建立文件，设文件名为 dzdt04.m 。 function dz=f(t,z) dz(1)=z(2); dz(2)=-2*t*z(4); dz(3)=z(4); dz(4)=2*t*z(2); dz=[dz(1);dz(2);dz(3);dz(4)]; （ 3 ）调用求解。 H=[-30,30];z0=[0;1;0;0]; [t,z]=ode45(dzdt04,H,z0); plot(z(:,1),z(:,3),m-) -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -2 -1.5 -1 -0.5 0 0.5 ? ? ? ? ? ? ? ? ? ? ? ? ? 0 2 0 2 x t y