科学工程计算与matlab编程5-3.pptVIP

5.2.4 微分方程转换5.2.4.1 单个高阶常微分方程处理方法 例： x0=[2;0]; t_final=3000; mu=1000; [t,y]=ode45(vdp_eq,[0,t_final],x0,[],mu); 由于变步长所采用的步长过小，所需时间较长，导致输出的y矩阵过大，超出计算机存储空间容量。所以不适合采用ode45()来求解，可用刚性方程求解算法ode15s( )。 6.2.4.2 高阶常微分方程组的变换方法 例： 描述函数： function dx=apolloeq(t,x) mu=1/82.45; mu1=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mu1)^2+x(3)^2); dx=[x(2); 2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3; x(4); -2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3]; 求解： x0=[1.2; 0; 0; -1;% 输入初值 tic, [t,y]=ode45(apolloeq,[0,20],x0); toc elapsed_time = 0.5000 length(t), plot(y(:,1),y(:,3)) ans = 689 得出的轨道不正确， 默认精度RelTol设置 得太大，从而导致的 误差传递，可减小该 值。 改变精度： options=odeset; options.RelTol=1e-6; tic, [t1,y1]=ode45(apolloeq,[0,20],x0,options); toc elapsed_time = 0.2970 length(t1), plot(y1(:,1),y1(:,3)), ans = 1873 min(diff(t1)) ans = 1.8927e-004 例： x0=[1.2; 0; 0; -1; tic, [t1,y1]=rk_4(apolloeq,[0,20,0.01],x0); toc elapsed_time = 0.8590 plot(y1(:,1),y1(:,3)) % 绘制出轨迹曲线 显而易见，这样求解 是错误的，应该采用 更小的步长。 tic, [t2,y2]=rk_4(apolloeq,[0,20,0.001],x0); toc elapsed_time = 12.4380 ％计算时间过长 plot(y2(:,1),y2(:,3)) % 绘制出轨迹曲线 严格说来某些点仍不 满足10－6的误差限， 所以求解常微分方程 组时建议采用变步长 算法，而不是定步长 算法。 例：试将二元方程组： 用MATLAB符号工具箱求解， 令 5.3特殊微分方程的数值解 5.3.1 刚性微分方程的求解 刚性微分方程 －刚性过程：微分方程描述的变化过程中，若包含着多 个相互作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有 “刚性”。 －刚性微分方程：描述“刚性”过程的微分方程称为刚性微分方程，相应的初值问题称为“刚性问题” 。 MATLAB采用求解函数ode15s()，该函数的调用格式和ode45()完全一致。 [t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,…) 例： ％计算 h_opt=odeset; h_opt.RelTol=1e-6; x0=[2;0]; t_final=3000; tic, mu=1000; [t,y]=ode15s(vdp_eq,[0,t_final],x0,h_opt,mu); toc elapsed_time = 1.2970 ％作图 plot(t,y(:,1)); figure; plot(t,y(:,2)) ％Vandpol方程输出结果y=(y,y’) y(:,1): 解曲线；部分曲线变化较平滑，某些点上变化在较快, y(:,2) ：解的变化率 例： 定义函数 function dy=c5exstf2(t,y) dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2; -10^4*y(1)+3000*(1-y(2))^2]; 方法一 t