第六章__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)] （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-) 例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 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-)