科学工程计算与matlab编程5-2

1、5.2 微分方程问题的数值解法5.2.1 算法概述0.( , ), ( )dxf t xatbdtx ax 微分方程求解的误差与步长问题：000121( )1,2,.,:nkkkkatx tattttbyknhtth时刻系统状态向量的值为数值解法：即求微分方程初值问题的解在若干点 的近似值 （）的方法。步长（一般取等步长 ）应用差商近似导数舍入误差10000( ,)xxhf txR0( )x t提高数值解精度的方法之一： 减小步长h的值； function outx,outy=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自变量的初值

2、和终值; %y0:函数在x0处的值,其可以为向量形式; %PointNum表示自变量在x0,xt上取的点数if nargin=4 | PointNumf1=(x,y)sin(x)+y; 2*pi/0.4 % 计算 xt=2pi 时应取得点数ans = 15.7080 x1,y1=MyEuler(f1,0,2*pi,1,16); %h=0.4 欧拉法所的的解x11,y11=MyEuler(f1,0,2*pi,1,32); % h=0.2 欧拉法所的的解h=0.2和h=0.4的数值解。y=dsolve(Dy=y+sin(t),y(0)=1); %该常微分方程的解析解for k=1:33t(k)=x

3、11(k);y2(k)=subs(y,t(k); %求其对应点的离散解endplot(x1,y1,+b,x11,y11,og,x11,y2,*r)legend(h=0.4的欧拉法解,h=0.2的欧拉解,符号解)误差判断：改进的Euler算法 ( , )dxf t xdt1 kkttdt11()( )kkkkx tx txx1 kkttdt11( ,)(,)2kkkkhf txf tx1111( ,)( ,)(,)2kkkkkkkkkkxxhf txhxxf txf tx改进Euler迭代公式Euler公式梯形公式function Xout,Yout=MyEulerPro(fun,x0,xt,y

4、0,PointNumber) %MyEulerPro 用改进的欧拉法解微分方程if nargin=4 | PointNumber5.2.2 四阶定步长Runge-Kutta算法及Matlab实现Euler算法：一阶Taylor多项式近似函数基本思想：2()O h误差：RK一般算法：推广：高级Taylor多项式近似函数提高数值解精度11111( ,)(,) (2,3,4,.)nniininijjpjnnjjkf txkf ta h xb kixxhkp为展开的项数例如：p4，即Matlab 实现：实现： function tout,yout=rk_4(odefile,tspan,y0) y0初值

5、列向量 t0=tspan(1); th=tspan(2); if length(tspan) comet3(x(:,1),x(:,2),x(:,3) 描述微分方程是常微分方程初值问题数值求解的关键。 f1=inline(-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);,. -x(1)*x(2)+28*x(2)-x(3),t,x); t_final=100; x0=0;0;1e-10; t,x=ode45(f1,0,t_final,x0); plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3); axis(10 42 -20 20

6、 -20 25);得出完全一致的结果。5.2.3.3 MATLAB 下带有附加参数的微分方程求解 例： 编写函数function xdot=lorenz1(t,x,flag,beta,rho,sigma) % flag变量是不能省略的 xdot=-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3); -x(1)*x(2)+sigma*x(2)-x(3);3/8,10,28 t_final=100; x0=0;0;1e-10; b2=3/8; r2=10; s2=28;%变量名不必一定为 等 t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2

7、,s2); plot(t2,x2), %options位置为，表示不需修改控制选项 figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(10 42 -20 20 -20 25);例：求解时的Lorenz微分方程 t_final=100; x0=0;0;1e-10; b2=2; r2=5; s2=20;%变量名不必一定为 等 t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2,s2); plot(t2,x2), %options位置为，表示不需修改控制选项 figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(0 72 -20 22 -35 40);f2=inline(-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3);,. -x(1)*x(2)+sigma*x(2)-x(3), t,x,flag,beta,rho,sigma); flag变量是不能省略t_final=100; x0=0;0;1e-10; b2=2; r2=5; s2=20;%变量名不必一定为 等 t2,x2=ode45(f2,0,t_final,x0,