matlab在数值分析中的应用 Rungekutta.pptVIP

求解： x0=[1.2; 0; 0; -1; tic, [t,y]=ode45(apolloeq,[0,20],x0); toc elapsed_time = 0.8310 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.8110 length(t1), plot(y1(:,1),y1(:,3)), ans = 1873 min(diff(t1)) ans = 1.8927e-004 plot(t1(1:end-1),… diff(t1)) 例： x0=[1.2; 0; 0; -1; tic, [t1,y1]=rk_4(apolloeq,[0,20,0.01],x0); toc elapsed_time = 4.2570 plot(y1(:,1),y1(:,3)) % 绘制出轨迹曲线 显而易见，这样求解 是错误的，应该采用 更小的步长。 tic, [t2,y2]=rk_4(apolloeq,[0,20,0.001],x0); toc elapsed_time = 124.4990 ％计算时间过长 plot(y2(:,1),y2(:,3)) % 绘制出轨迹曲线 严格说来某些点仍不 满足10－6的误差限， 所以求解常微分方程 组时建议采用变步长 算法，而不是定步长 算法。 例： 用MATLAB符号工具箱求解， 令 ％ syms x1 x2 x3 x4 [dx,dy]=solve(dx+2*x4*x1=2*dy, dx*x4+ … 3*x2*dy+x1*x4-x3=5,dx,dy) dx = -2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2) dy = (2*x4^2*x1-x4*x1+x3+5)/(2*x4+3*x2) 对于更复杂的问题来说，手工变换的难度将很大，所以如有可能，可采用计算机去求解有关方程，获得解析解。如不能得到解析解，也需要在描写一阶常微分方程组时列写出式子，得出问题的数值解。 7.3特殊微分方程的数值解 7.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 = 2.5240 ％作图 plot(t,y(:,1)); figure; plot(t,y(:,2)) y(:,1)曲线变化较平滑, y(:,2)变化在某些点上较快。 例： 定义函数 function dy=c7exstf2(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]; 方法一 tic,[t2,y2]=ode45(c7exstf2,[0,100],[0;1]); toc elapsed_time = 229.4700 length(t2), plot(t2,y2) ans = 356941 步长分析： format long, [min(diff(t2)), max(diff(t2))] ans = 0.00022220693884 0.00214971787184 plot(t2(1:end-1),diff(t2)) 方法二，用ode15s()代替ode