matlab在数值分析中的应用7 .ppt

求解： 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的误差限， 所以求解常微分方程 组时建议采用变步长 算法，而不是定步长 算法。 第七章 微分方程问题的解法 微分方程的解析解方法 微分方程问题的数值解法 微分方程问题算法概述 四阶定步长 Runge-Kutta算法及 MATLAB 实现 一阶微分方程组的数值解 微分方程转换 特殊微分方程的数值解 7.1 微分方程的解析解方法 格式： y=dsolve(f1, f2, …, fm) 格式：指明自变量 y=dsolve(f1, f2, …, fm ,’x’) fi即可以描述微分方程，又可描述初始条件或边界条件。如： 描述微分方程时 描述条件时 例： syms t; u=exp(-5*t)*cos(2*t+1)+5; uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu = 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10 syms t y; y=dsolve([D4y+10*D3y+35*D2y+50*Dy+24*y=,... 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10]) y=dsolve([D4y+10*D3y+35*D2y+50*Dy+24*y=,... 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) ... +10], y(0)=3, Dy(0)=2, D2y(0)=0, D3y(0)=0) 分别处理系数，如： [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))] ans = -8704 185 % rat()最接近有理数的分数 判断误差： vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185) ans = .114731975864790922564144636e-4 y=dsolve([D4y+10*D3y+35*D2y+50*Dy+24*y=,... 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) + ... 10],y(0)=1/2,Dy(pi)=1,D2y(2*pi)=0,Dy(2*pi)=1/5); 如果用推导的方法求Ci的值，每个系数的解析解至少要写出10数行,故可采用有理式近似 的方式表示. vpa(y,10) %有理近似值 ans = 1.196361839*exp(-5.*t