matlab在数值分析中的应用

第七章微分方程问题的解法•微分方程的解析解方法•微分方程问题的数值解法–微分方程问题算法概述–四阶定步长Runge-Kutta算法及MATLAB实现–一阶微分方程组的数值解–微分方程转换•特殊微分方程的数值解7.1微分方程的解析解方法•格式：y=dsolve(f1,f2,…,fm)•格式：指明自变量y=dsolve(f1,f2,…,fm,’x’)fi即可以描述微分方程，又可描述初始条件或边界条件。如：描述微分方程时描述条件时(4)()747ytDy(2)32(2)3yDy例：symst;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*uuu=87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10symsty;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=-8704185%rat()最接近有理数的分数判断误差：vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)ans=.114731975864790922564144636e-4y=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)+.4166666667-.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e-1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(-5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(-3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(-1.*t)+442590.9059*exp(-4.*t)•例：symstxx=dsolve('Dx=x*(1-x^2)')x=[1/(1+exp(-2*t)*C1)^(1/2)][-1/(1+exp(-2*t)*C1)^(1/2)]symstx;x=dsolve('Dx=x*(1-x^2)+1')Warning:Explicitsolutioncouldnotbefound;implicitsolutionreturned.InD:\MATLAB6p5\toolbox\symbolic\dsolve.matline292x=t-Int(1/(a-a^3+1),a=``..x)+C1=0故只有部分非线性微分方程有解析解。7.2微分方程问题的数值解法7.2.1微分方程问题算法概述微分方程求解的误差与步长问题：7.2.2四阶定步长Runge-Kutta算法及MATLAB实现function[tout,yout]=rk_4(odefile,tspan,y0)％y0初值列向量t0=tspan(1);th=tspan(2);iflength(tspan)=3,h=tspan(3);%tspan=[t0,th,h]else,h=tspan(2)-tspan(1);th=tspan(end);end％等间距数组tout=[t0:h:th]';yout=[];fort=tout'k1=h*eval([odefile‘(t,y0)’]);%odefile是一个字符串变量，为表示微分方程的文件名。k2=h*eval([odefile'(t+h/2,y0+0.5*k1)']);k3=h*eval([odefile'(t+h/2,y0+0.5*k2)']);k4=h*eval([odefile'(t+h,y0+k3)']);y0=y0+(k1+2*k2+2*k3+k4)/6;yout=[yout;y0'];end％由效果看，该算法不是一个较好的方法。7.2.3一阶微分方程组的数值解7.2.3.1四阶五级Runge-Kutta-Felhberg算法通过误差向量调节步长，此为自动变步长方法。四阶五级RKF算法有参量系数表。7.2.3.2基于MATLAB的微分方程求解函数格式1：直接求解[t,x]=ode45(Fun,[t0,tf],x0)格式2：带有控制参数[t,x]=ode45(Fun,[t0,tf],x0，options)格式3：带有附加参数[t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,…)[t0,tf]求解区间，x0初值问题的初始状态变量。描述需要求解的微分方程组：不需附加变量的格式functionxd=funname(t,x)可以使用附加变量functionxd=funname(t,x,flag,p1,p2,…)％t是时间变量或自变量（必须给），x为状态向量，xd为状态向量的导数。flag用来控制求解过程,指定初值，即使初值不用指定，也必须有该变量占位。修改变量：options唯一结构体变量，用odeset()修改。options=odeset(‘RelTol’,1e-7);options=odeset;options.RelTol=1e-7;•例：自变函数functionxdot=lorenzeq(t,x)xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);…-x(1)*x(2)+28*x(2)-x(3)];t_final=100;x0=[0;0;1e-10];%t_final为设定的仿真终止时间[t,x]=ode45('lorenzeq',[0,t_final],x0);plot(t,x),figure;%打开新图形窗口plot3(x(:,1),x(:,2),x(:,3));axis([1042-2020-2025]);%根据实际数值手动设置坐标系•可采用comet3()函数绘制动画式的轨迹。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([1042-2020-2025]);得出完全一致的结果。7.2.3.3MATLAB下带有附加参数的微分方程求解•例：•编写函数functionxdot=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)];求微分方程：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([072-2022-3540]);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变量是不能省略的7.2.4微分方程转换7.2.4.1单个高阶常微分方程处理方法•例：•函数描述为：functiony=vdp_eq(t,x,flag,mu)y=[x(2);-mu*(x(1)^2-1)*x(2)-x(1)];x0=[-0.2;-0.7];t_final=20;mu=1;[t1,y1]=ode45('vdp_eq',[0,t_final],x0,[],mu);mu=2;[t2,y2]=ode45('vdp_eq',[0,t_final],x0,[],mu);plot(t1,y1,t2,y2,':')figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),':')x0=[2;0];t_final=3000;mu=1000;[t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu);由于变步长所采用的步长过小，所需时间较长，导致输出的y矩阵过大，超出计算机存储空间容量。所以不适合采用ode45()来求解，可用刚性方程求解算法ode15s()。7.2.4.2高阶常微分方程组的变换方法•例：•描述函数：functiondx=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.04935751];tic,[t,y]=ode45('apolloeq',[0,20],x0);tocelapsed_time=0.8310length(t),plot(y(:,1),y(:,3))ans=689得出的轨道不正确，默认精度RelTol设置得太大，从而导致的误差传递，可减小该值。•改变精度：options=odeset;options.RelTol=1e-6;tic,[t1,y1]=ode45('apolloeq',[0,20],x0,options);tocelapsed_time=0.8110length(t1),plot(y1(:,1),y1(:,3)),ans=1873min(diff(t1))ans=1.8927e-004plot(t1(1:end-1),…diff(t1))•例：x0=[1.2;0;0;-1.04935751];tic,[t1,y1]=rk_4('apolloeq',[0,20,0.01],x0);tocelapsed_time=4.2570plot(y1(:,1),y1(:,3))%绘制出轨迹曲线显而易见，这样求解是错误的，应该采用更小的步长。tic,