微分方程问题的解法学习教案.pptx

会计学1微分方程微分方程(wi fn fn chn)问题的解问题的解法法第一页，共85页。y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*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).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)+10,y(0)=3,Dy(0)=2,D2y(0)=0,D3y(0)=0)第2页/共85页第二页，共85页。分别处理系数，如：分别处理系数，如：n,d=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)n,d=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)ans=ans=-8704 185%rat()-8704 185%rat()最接近最接近(jijn)(jijn)有理数的分数有理数的分数判断误差：判断误差：vpa(-445/26*cos(sym(1)-51/13*sin(1)-69/2+8704/185)vpa(-445/26*cos(sym(1)-51/13*sin(1)-69/2+8704/185)ans=ans=第3页/共85页第三页，共85页。y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*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)+.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);10,y(0)=1/2,Dy(pi)=1,D2y(2*pi)=0,Dy(2*pi)=1/5);如果用推导的方法求如果用推导的方法求CiCi的值，每个系数的解析解至少要的值，每个系数的解析解至少要写出写出1010数行数行,故可采用有理式近似故可采用有理式近似 的方式表示的方式表示(bi(bi osh).osh).vpa(y,10)%vpa(y,10)%有理近似值有理近似值ans=ans=1.196361839*exp(-5.*t)+.4166666667-.4785447354*sin(t)*cos(t)*exp1.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)-.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(-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)(-2.*t)-219.1293619*exp(-1.*t)+442590.9059*exp(-4.*t)第4页/共85页第四页，共85页。n n例：求解(qi ji)n n x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t),Dy=4*x+3*y+4*exp(-t)第5页/共85页第五页，共85页。n n例：例：n n syms t x syms t xn n x=dsolve(Dx=x*(1-x2)x=dsolve(Dx=x*(1-x2)n nx=x=n n 1/(1+exp(-2*t)*C1)(1/2)1/(1+exp(-2*t)*C1)(1/2)n n-1/(1+exp(-2*t)*C1)(1/2)-1/(1+exp(-2*t)*C1)(1/2)n n syms t x;x=dsolve(Dx=x*(1-x2)+1)syms t x;x=dsolve(Dx=x*(1-x2)+1)n nWarning:Explicit solution could not be found;implicit Warning:Explicit solution could not be found;implicit solution returned.solution returned.n n In D:MATLAB6p5toolboxsymbolicdsolve.m at In D:MATLAB6p5toolboxsymbolicdsolve.m at line 292line 292n nx=x=n nt-Int(1/(a-a3+1),a=.x)+C1=0t-Int(1/(a-a3+1),a=.x)+C1=0n n故只有故只有(zh(zh y y u)u)部分非线性微分方程有解析解。部分非线性微分方程有解析解。第6页/共85页第六页，共85页。6.2 6.2 微分方程问题的数值解法微分方程问题的数值解法微分方程问题算法微分方程问题算法(sun f)(sun f)概述概述第7页/共85页第七页，共85页。微分方程(wi fn fn chn)求解的误差与步长问题：第8页/共85页第八页，共85页。第9页/共85页第九页，共85页。第10页/共85页第十页，共85页。四阶定步长四阶定步长Runge-KuttaRunge-Kutta算法算法(sun f)(sun f)及及 MATLAB MATLAB 实现实现第11页/共85页第十一页，共85页。function tout,yout=rk_4(odefile,tspan,y0)function tout,yout=rk_4(odefile,tspan,y0)y0y0初值列向量初值列向量(xingling)(xingling)t0=tspan(1);th=tspan(2);t0=tspan(1);th=tspan(2);if length(tspan)=3,h=tspan(3);%tspan=t0,th,h if length(tspan)t_final=100;x0=0;0;1e-10;%t_final t_final=100;x0=0;0;1e-10;%t_final为设定的仿真为设定的仿真(f(f n n zhn)zhn)终止时间终止时间 t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure;%figure;%打开新图形窗口打开新图形窗口 plot3(x(:,1),x(:,2),x(:,3);plot3(x(:,1),x(:,2),x(:,3);axis(10 42-20 20-20 25);%axis(10 42-20 20-20 25);%根据实际数值手动设置坐标根据实际数值手动设置坐标系系第17页/共85页第十七页，共85页。n n可采用(ciyng)comet3()函数绘制动画式的轨迹。n n comet3(x(:,1),x(:,2),x(:,3)第18页/共85页第十八页，共85页。n n描述微分方程是常微分方程初值问题数值求解的关键(gunjin)。n n f1=inline(-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);,.n n -x(1)*x(2)+28*x(2)-x(3),t,x);n n t_final=100;x0=0;0;1e-10;n n t,x=ode45(f1,0,t_final,x0);n n plot(t,x),figure;n n plot3(x(:,1),x(:,2),x(:,3);axis(10 42-20 20-20 25);n n得出完全一致的结果。第19页/共85页第十九页，共85页。下带有附加参数的微分方程下带有附加参数的微分方程下带有附加参数的微分方程下带有附加参数的微分方程(wi fn(wi fn(wi fn(wi fn fn chn)fn chn)fn chn)fn chn)求解求解求解求解n n例：第20页/共85页第二十页，共85页。n n编写函数编写函数n nfunction xdot=lorenz1(t,x,flag,beta,rho,sigma)function xdot=lorenz1(t,x,flag,beta,rho,sigma)n n flag flag变量变量(binling)(binling)是是不能省略的不能省略的n n xdot=-beta*x(1)+x(2)*x(3);xdot=-beta*x(1)+x(2)*x(3);n n -rho*x(2)+rho*x(3);-rho*x(2)+rho*x(3);n n -x(1)*x(2)+sigma*x(2)-x(3);-x(1)*x(2)+sigma*x(2)-x(3);n n求微分方程：求微分方程：n n t_final=100;x0=0;0;1e-10;t_final=100;x0=0;0;1e-10;n n b2=2;r2=5;s2=20;b2=2;r2=5;s2=20;n n t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2,s2);t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2,s2);n n plot(t2,x2),plot(t2,x2),optionsoptions位置为位置为，表示不需修改控制选项，表示不需修改控制选项n n figure;plot3(x2(:,1),x2(:,2),x2(:,3);axis(0 72-20 22 -35 40);figure;plot3(x2(:,1),x2(:,2),x2(:,3);axis(0 72-20 22 -35 40);第21页/共85页第二十一页，共85页。n nf2=inline(-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);,.n n -x(1)*x(2)+sigma*x(2)-x(3),t,x,flag,beta,rho,sigma);n n flag变量(binling)是不能省略的第22页/共85页第二十二页，共85页。微分方程微分方程(wi fn fn chn)(wi fn fn chn)转换转换单个高阶常微分方程单个高阶常微分方程(wi fn fn(wi fn fn chn)chn)处理方法处理方法第23页/共85页第二十三页，共85页。第24页/共85页第二十四页，共85页。n n例：例：n n函数函数(hnsh)(hnsh)描述为：描述为：n n function y=vdp_eq(t,x,flag,mu)function y=vdp_eq(t,x,flag,mu)n n y=x(2);-mu*(x(1)2-1)*x(2)-x(1);y=x(2);-mu*(x(1)2-1)*x(2)-x(1);n n x0=-0.2;-0.7;t_final=20;x0=-0.2;-0.7;t_final=20;n n mu=1;t1,y1=ode45(vdp_eq,0,t_final,x0,mu);mu=1;t1,y1=ode45(vdp_eq,0,t_final,x0,mu);n n mu=2;t2,y2=ode45(vdp_eq,0,t_final,x0,mu);mu=2;t2,y2=ode45(vdp_eq,0,t_final,x0,mu);n n plot(t1,y1,t2,y2,:)plot(t1,y1,t2,y2,:)n n figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:)figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:)第25页/共85页第二十五页，共85页。n n x0=2;0;t_final=3000;x0=2;0;t_final=3000;n n mu=1000;t,y=ode45(vdp_eq,0,t_final,x0,mu);mu=1000;t,y=ode45(vdp_eq,0,t_final,x0,mu);n n 由于变步长所采用的步长过小，所需时间较长，由于变步长所采用的步长过小，所需时间较长，导致输出的导致输出的y y矩阵过大，超出计算机存储空间容量。矩阵过大，超出计算机存储空间容量。所以不适合采用所以不适合采用ode45()ode45()来求解，可用刚性方程来求解，可用刚性方程(fngchng)(fngchng)求解算法求解算法ode15s()ode15s()。第26页/共85页第二十六页，共85页。高阶常微分方程组的变换高阶常微分方程组的变换(binhun)(binhun)方法方法第27页/共85页第二十七页，共85页。n n例：第28页/共85页第二十八页，共85页。第29页/共85页第二十九页，共85页。n n描述(mio sh)函数：n n function dx=apolloeq(t,x)n n mu=1/82.45;mu1=1-mu;n n r1=sqrt(x(1)+mu)2+x(3)2);n n r2=sqrt(x(1)-mu1)2+x(3)2);n n dx=x(2);n n 2*x(4)+x(1)-mu1*(x(1)+mu)/r13-mu*(x(1)-mu1)/r23;n n x(4);n n -2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;第30页/共85页第三十页，共85页。n n求解：求解：n n x0=1.2;0;0;-1.04935751;x0=1.2;0;0;-1.04935751;n n tic,t,y=ode45(apolloeq,0,20,x0);toc tic,t,y=ode45(apolloeq,0,20,x0);tocn nelapsed_time=elapsed_time=n n 0.8310 0.8310n n length(t),length(t),n n plot(y(:,1),y(:,3)plot(y(:,1),y(:,3)n nans=ans=n n 689 689n n得出的轨道不正确得出的轨道不正确(zhngqu)(zhngqu)，n n默认精度默认精度RelTolRelTol设置设置n n得太大，从而导致的得太大，从而导致的n n误差传递，可减小该误差传递，可减小该n n值。值。第31页/共85页第三十一页，共85页。n n改变(gibin)精度：n n options=odeset;options.RelTol=1e-6;n n tic,t1,y1=ode45(apolloeq,0,20,x0,options);tocn nelapsed_time=n n 0.8110n n length(t1),n n plot(y1(:,1),y1(:,3),n nans=n n 1873第32页/共85页第三十二页，共85页。min(diff(t1)min(diff(t1)ans=ans=1.8927e-004 1.8927e-004 plot(t1(1:end-1),plot(t1(1:end-1),diff(t1)diff(t1)第33页/共85页第三十三页，共85页。n n例：第34页/共85页第三十四页，共85页。n n x0=1.2;0;0;-1.04935751;n n tic,t1,y1=rk_4(apolloeq,0,20,0.01,x0);tocn nelapsed_time=n n 4.2570n n plot(y1(:,1),y1(:,3)n n%绘制出轨迹曲线n n显而易见，这样(zhyng)求解n n是错误的，应该采用n n更小的步长。第35页/共85页第三十五页，共85页。n n tic,t2,y2=rk_4(apolloeq,0,20,0.001,x0);toc tic,t2,y2=rk_4(apolloeq,0,20,0.001,x0);tocn nelapsed_time=elapsed_time=n n 124.4990 124.4990 计算时间过长计算时间过长n n plot(y2(:,1),y2(:,3)plot(y2(:,1),y2(:,3)n n%绘制出轨迹曲线绘制出轨迹曲线n n严格说来某些点仍不严格说来某些点仍不n n满足满足10106 6的误差的误差(wch)(wch)限，限，n n所以求解常微分方程所以求解常微分方程n n组时建议采用变步长组时建议采用变步长n n算法，而不是定步长算法，而不是定步长n n算法。算法。第36页/共85页第三十六页，共85页。n n例：第37页/共85页第三十七页，共85页。用用MATLABMATLAB符号工具箱求解，符号工具箱求解，令令 syms x1 x2 x3 x4 syms x1 x2 x3 x4 dx,dy=solve(dx+2*x4*x1=2*dy,dx*x4+dx,dy=solve(dx+2*x4*x1=2*dy,dx*x4+3*x2*dy+x1*x4-x3=5,dx,dy)%dx,dy3*x2*dy+x1*x4-x3=5,dx,dy)%dx,dy为指定变量为指定变量dx=dx=-2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2)-2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2)dy=dy=(2*x42*x1-x4*x1+x3+5)/(2*x4+3*x2)(2*x42*x1-x4*x1+x3+5)/(2*x4+3*x2)对于更复杂的问题来说，手工变换的难度将很大，所对于更复杂的问题来说，手工变换的难度将很大，所以如有可能，可采用计算机去求解有关方程，获得解析以如有可能，可采用计算机去求解有关方程，获得解析解。如不能得到解析解，也需要解。如不能得到解析解，也需要(xyo)(xyo)在描写一阶常微在描写一阶常微分方程组时列写出式子，得出问题的数值解。分方程组时列写出式子，得出问题的数值解。第38页/共85页第三十八页，共85页。6.36.3特殊微分方程的数值特殊微分方程的数值(shz)(shz)解解 刚性微分方程的求解刚性微分方程的求解n n刚性微分方程n n 一类特殊的常微分方程，其中一些解变化缓慢，另一些变化快，且相差悬殊，这类方程常常称为刚性方程。n nMATLAB采用求解函数(hnsh)ode15s()，该函数(hnsh)的调用格式和ode45()完全一致。n nt,x=ode15s(Fun,t0,tf,x0,options,p1,p2,)第39页/共85页第三十九页，共85页。n n例：n n计算(j sun)n n h_opt=odeset;h_opt.RelTol=1e-6;n n x0=2;0;t_final=3000;n n tic,mu=1000;t,y=ode15s(vdp_eq,0,t_final,x0,h_opt,mu);tocn nelapsed_time=n n 2.5240第40页/共85页第四十页，共85页。作图 plot(t,y(:,1);figure;plot(t,y(:,2)y(:,1)曲线(qxin)变化较平滑,y(:,2)变化在某些点上较快。第41页/共85页第四十一页，共85页。n n例：n n定义(dngy)函数n nfunction dy=c7exstf2(t,y)n n dy=0.04*(1-y(1)-(1-y(2)*y(1)+0.0001*(1-y(2)2;n n -104*y(1)+3000*(1-y(2)2;第42页/共85页第四十二页，共85页。n n方法(fngf)一n n tic,t2,y2=ode45(c7exstf2,0,100,0;1);tocn nelapsed_time=n n 229.4700n n length(t2),plot(t2,y2)n nans=n n 356941第43页/共85页第四十三页，共85页。n n步长分析(fnx)：n n format long,min(diff(t2),max(diff(t2)n nans=n n plot(t2(1:end-1),diff(t2)第44页/共85页第四十四页，共85页。n n方法(fngf)二，用ode15s()代替ode45()n n opt=odeset;opt.RelTol=1e-6;n n tic,t1,y1=ode15s(c7exstf2,0,100,0;1,opt);tocn nelapsed_time=n n length(t1),n n plot(t1,y1)n nans=n n 169第45页/共85页第四十五页，共85页。隐式微分方程隐式微分方程隐式微分方程隐式微分方程(wi fn fn chn)(wi fn fn chn)(wi fn fn chn)(wi fn fn chn)求解求解求解求解n n隐式微分方程隐式微分方程(wi fn fn(wi fn fn chn chn)为不能转化为显式常为不能转化为显式常微分方程微分方程(wi fn fn(wi fn fn chn chn)组的方程组的方程n n例：例：第46页/共85页第四十六页，共85页。n n编写函数编写函数(hnsh)(hnsh)：n nfunction dx=c7ximp(t,x)function dx=c7ximp(t,x)n n A=sin(x(1)cos(x(2);-cos(x(2)sin(x(1);A=sin(x(1)cos(x(2);-cos(x(2)sin(x(1);n n B=1-x(1);-x(2);dx=inv(A)*B;B=1-x(1);-x(2);dx=inv(A)*B;n n求解：求解：n n opt=odeset;opt.RelTol=1e-6;opt=odeset;opt.RelTol=1e-6;n n t,x=ode45(c7ximp,0,10,0;0,opt);plot(t,x)t,x=ode45(c7ximp,0,10,0;0,opt);plot(t,x)第47页/共85页第四十七页，共85页。微分代数方程微分代数方程(dish fngchng)(dish fngchng)求求解解例：第48页/共85页第四十八页，共85页。n n编写函数n n function dx=c7eqdae(t,x)n n dx=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);n n 2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);n n x(1)+x(2)+x(3)-1;n n M=1,0,0;0,1,0;0,0,0;n n options=odeset;n n options.Mass=M;n n Mass微分代数方程中n n的质量(zhling)矩阵（控制参数）n n x0=0.8;0.1;0.1;n n t,x=ode15s(c7eqdae,0,20,x0,options);plot(t,x)第49页/共85页第四十九页，共85页。编写(binxi)函数：function dx=c7eqdae1(t,x)dx=-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2);第50页/共85页第五十页，共85页。x0=0.8;0.1;x0=0.8;0.1;fDae=inline(-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);,.fDae=inline(-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);,.2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2),t,x);2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2),t,x);t1,x1=ode45(fDae,0,20,x0);plot(t1,x1,t1,1-sum(x1)t1,x1=ode45(fDae,0,20,x0);plot(t1,x1,t1,1-sum(x1)第51页/共85页第五十一页，共85页。延迟微分方程延迟微分方程(wi fn fn chn)(wi fn fn chn)求解求解sol:结构体数据，sol.x:时间(shjin)向量t,sol.y:状态向量。第52页/共85页第五十二页，共85页。n n例：例：第53页/共85页第五十三页，共85页。编写函数编写函数(hnsh)(hnsh)：function dx=c7exdde(t,x,z)function dx=c7exdde(t,x,z)xlag1=z(:,1);%xlag1=z(:,1);%第一列表示提取第一列表示提取 xlag2=z(:,2);xlag2=z(:,2);dx=1-3*x(1)-xlag1(2)-0.2*xlag2(1)3-xlag2(1);dx=1-3*x(1)-xlag1(2)-0.2*xlag2(1)3-xlag2(1);x(3);4*x(1)-2*x(2)-3*x(3);x(3);4*x(1)-2*x(2)-3*x(3);历史数据函数历史数据函数(hnsh)(hnsh)：function S=c7exhist(t)function S=c7exhist(t)S=zeros(3,1);S=zeros(3,1);第54页/共85页第五十四页，共85页。n n求解：求解：n n lags=1 0.5;tx=dde23(c7exdde,lags,zeros(3,1),0,10);lags=1 0.5;tx=dde23(c7exdde,lags,zeros(3,1),0,10);n n plot(tx.x,tx.y(2,:)plot(tx.x,tx.y(2,:)n n与与ode45()ode45()等返回的等返回的x x矩阵矩阵(j(j zhn)zhn)不一样，它是按行不一样，它是按行排列的。排列的。第55页/共85页第五十五页，共85页。6.46.4边值问题的计算机求解边值问题的计算机求解(qi(qi ji)ji)第56页/共85页第五十六页，共85页。边值问题的打靶边值问题的打靶(d b)(d b)算法算法数学方法描述数学方法描述(mio sh)(mio sh)：以二阶方程为例以二阶方程为例 第57页/共85页第五十七页，共85页。编写编写(binxi)(binxi)函数：函数：线性的线性的function t,y=shooting(f1,f2,tspan,x0f,varargin)function t,y=shooting(f1,f2,tspan,x0f,varargin)t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);t,y1=ode45(f1,tspan,1;0,varargin);t,y1=ode45(f1,tspan,1;0,varargin);t,y2=ode45(f1,tspan,0;1,varargin);t,y2=ode45(f1,tspan,0;1,varargin);t,yp=ode45(f2,tspan,0;0,varargin);t,yp=ode45(f2,tspan,0;0,varargin);m=(gb-ga*y1(end,1)-yp(end,1)/y2(end,1);m=(gb-ga*y1(end,1)-yp(end,1)/y2(end,1);t,y=ode45(f2,tspan,ga;m,varargin);t,y=ode45(f2,tspan,ga;m,varargin);第58页/共85页第五十八页，共85页。例：例：编写编写(binxi)(binxi)函数：函数：function xdot=c7fun1(t,x)function xdot=c7fun1(t,x)xdot=x(2);-2*x(1)+3*x(2);xdot=x(2);-2*x(1)+3*x(2);function xdot=c7fun2(t,x)function xdot=c7fun2(t,x)xdot=x(2);t-2*x(1)+3*x(2);xdot=x(2);t-2*x(1)+3*x(2);t,y=shooting(c7fun1,t,y=shooting(c7fun1,c7fun2,0,1,1;2);plot(t,y)c7fun2,0,1,1;2);plot(t,y)第59页/共85页第五十九页，共85页。原方程原方程(fngchng)(fngchng)的解析解为的解析解为解的检验解的检验 y0=(exp(2)-3)*exp(t)+(3-exp(1)*exp(2*t)/(4*exp(1)*(exp(1)-y0=(exp(2)-3)*exp(t)+(3-exp(1)*exp(2*t)/(4*exp(1)*(exp(1)-1)+3/4+t/2;1)+3/4+t/2;norm(y(:,1)-y0)%norm(y(:,1)-y0)%整个解函数检验整个解函数检验ans=ans=4.4790e-008 4.4790e-008 norm(y(end,1)-2)%norm(y(end,1)-2)%终点条件检验终点条件检验ans=ans=2.2620e-008 2.2620e-008第60页/共85页第六十页，共85页。非线性方程非线性方程(fngchng)(fngchng)边值问题的打靶算法：边值问题的打靶算法：用用NewtonNewton迭代法处理迭代法处理第61页/共85页第六十一页，共85页。编写编写(binxi)(binxi)函数：函数：function t,y=nlbound(funcn,funcv,tspan,x0f,tol,varargin)function t,y=nlbound(funcn,funcv,tspan,x0f,tol,varargin)t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);m=1;t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);m=1;m0=0;m0=0;while(norm(m-m0)tol),m0=m;while(norm(m-m0)tol),m0=m;t,v=ode45(funcv,tspan,ga;m;0;1,varargin);t,v=ode45(funcv,tspan,ga;m;0;1,varargin);m=m0-(v(end,1)-gb)/(v(end,3);m=m0-(v(end,1)-gb)/(v(end,3);end end t,y=ode45(funcn,tspan,ga;m,varargin);t,y=ode45(funcn,tspan,ga;m,varargin);第62页/共85页第六十二页，共85页。例：例：编写两个编写两个(li(li n n )函数：函数：function xdot=c7fun3(t,x)function xdot=c7fun3(t,x)xdot=x(2);2*x(1)*x(2);x(4);2*x(2)*x(3)+2*x(1)*x(4);xdot=x(2);2*x(1)*x(2);x(4);2*x(2)*x(3)+2*x(1)*x(4);function xdot=c7fun4(t,x)function xdot=c7fun4(t,x)xdot=x(2);2*x(1)*x(2);xdot=x(2);2*x(1)*x(2);第63页/共85页第六十三页，共85页。t,y=nlbound(c7fun4,c7fun3,0,pi/2,-1,1,1e-8);t,y=nlbound(c7fun4,c7fun3,0,pi/2,-1,1,1e-8);plot(t,y);set(gca,xlim,0,pi/2);plot(t,y);set(gca,xlim,0,pi/2);精确精确(jngqu)(jngqu)解：解：检验：检验：y0=tan(t-pi/4);y0=tan(t-pi/4);norm(y(:,1)-y0)norm(y(:,1)-y0)ans=ans=1.6629e-005 1.6629e-005 norm(y(end,1)-1)norm(y(end,1)-1)ans=ans=5.2815e-006 5.2815e-006第64页/共85页第六十四页，共85页。线性微分方程的有限线性微分方程的有限线性微分方程的有限线性微分方程的有限(yuxin)(yuxin)(yuxin)(yuxin)差分算法差分算法差分算法差分算法把等式把等式(dngsh)(dngsh)左边用差商表示。左边用差商表示。第65页/共85页第六十五页，共85页。第66页/共85页第六十六页，共85页。编写编写(binxi)(binxi)函数：函数：function x,y=fdiff(funcs,tspan,x0f,n)function x,y=fdiff(funcs,tspan,x0f,n)t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);h=(tfinal-t0)/n;h=(tfinal-t0)/n;for i=1:n,x(i)=t0+h*(i-1);end,for i=1:n,x(i)=t0+h*(i-1);end,x0=x(1:n-1);t=-2+h2*feval(funcs,x0,2);x0=x(1:n-1);t=-2+h2*feval(funcs,x0,2);tmp=feval(funcs,x0,1);tmp=feval(funcs,x0,1);v=1+h*tmp/2;w=1-h*tmp/2;b=h2*feval(funcs,x0,3);v=1+h*tmp/2;w=1-h*tmp/2;b=h2*feval(funcs,x0,3);b(1)=b(1)-w(1)*ga;b(n-1)=b(n-1)-v(n-1)*gb;b=b;b(1)=b(1)-w(1)*ga;b(n-1)=b(n-1)-v(n-1)*gb;b=b;A=diag(t);A=diag(t);for i=1:n-2,A(i,i+1)=v(i);A(i+1,i)=w(i+1);end for i=1:n-2,A(i,i+1)=v(i);A(i+1,i)=w(i+1);end y=inv(A)*b;x=x tfinal;y=ga;y;gb;y=inv(A)*b;x=x tfinal;y=ga;y;gb;第67页/共85页第六十七页，共85页。例：编写(binxi)函数：function y=c7fun5(x,key)switch key case 1,y=1+x;case 2,y=1-x;otherwise,y=1+x.2;end t,y=fdiff(c7fun5,0,1,1,4,50);plot(t,y)第68页/共85页第六十八页，共8