数值分析常微分方程初值问题精品文稿.ppt
数值分析常微分方程初值问题第1页，本讲稿共17页例1.罗伦茨模型求解区域求解区域:0 t 80初始条件：初始条件：x(0)=0,y(0)=0,z(0)=0.0001取取 =8/3,=28,=102/16第2页，本讲稿共17页例2.Lotka-Volterra捕食者与被捕食者模型取取 a=1,b=0.01,c=1,d=0.02求解区域求解区域:0 t 15初始条件初始条件:x(0)=20,y(0)=203/16-x(t)-y(t)第3页，本讲稿共17页例例3.人口模型人口模型(Logistic):P=rP(1 P/K)取取:r=0.02,K=500初始条件初始条件:P(0)=76 例例4.传染病模型传染病模型:P=rP(L P)L=25000,r=0.00001 初始条件初始条件:P(0)=24/16第4页，本讲稿共17页常微分方程初值问题:例 数值方法：数值方法：0=x0 x1 x2 xN 2求未知函数求未知函数 y=y(x)x(x0,)求求:y(x1),y(x2),y(x3),y(xN)的近似值的近似值 y1,y2,y3,yN5/16第5页，本讲稿共17页连续问题:离散问题:n=1,2,3,6/16第6页，本讲稿共17页求解常微分方程初值问题的Euler方法 7/16第7页，本讲稿共17页例 取步长取步长 h=2/10,2/20,2/30,2/40,用欧拉法求解用欧拉法求解.解解:记记 f(x,y)=y x y2,xn=0.2n (n=0,1,2,N)由由Euler公式公式:yn+1=yn+h(yn xn yn2)(n=0,1,N)y0=1 N 10 20 30 40 h 0.2 0.1 0.0667 0.05误差 0.1059 0.0521 0.0342 0.02568/16第8页，本讲稿共17页n=input(input n:=);f=inline(y-x.*y.2);h=2/n;x=h:h:2;y0=1;y(1)=y0+h*f(0,y0);for k=1:n-1 xk=x(k);yk=y(k);y(k+1)=yk+h*f(xk,yk);endt=-1:h/2:5;yy=1./(t-1+2*exp(-t);plot(0,1,ro,x,y,ro,t,yy)y1=1./(x-1+2*exp(-x);error=max(abs(y-y1)o:数值解数值解-:准确解准确解 9/16解析解:第9页，本讲稿共17页y=f(x,y)梯形公式梯形公式:左矩形公式左矩形公式用数值积分方法离散化常微分方程10/16第10页，本讲稿共17页算法算法:k1=f(xn,yn),k2=f(xn+1,yn+h k1)由梯形公式推出的预-校方法:11/16第11页，本讲稿共17页n=input(input n:=);f=inline(y-x.*y.2);h=2/n;x=h:h:2;y0=1;k1=f(0,y0);k2=f(h,y0+h*k1);y(1)=y0+0.5*h*(k1+k2);for k=1:n-1 xk=x(k);yk=y(k);k1=f(xk,yk);k2=f(xk+h,yk+h*k1);y(k+1)=yk+0.5*h*(k1+k2);endt=-1:h/2:5;yy=1./(t-1+2*exp(-t);plot(0,1,ro,x,y,ro,t,yy)y1=1./(x-1+2*exp(-x);error=max(abs(y-y1)input n:=10error=0.012312/16 o:数值解数值解-:准确解准确解 第12页，本讲稿共17页预预-校方法校方法,h=0.2时时误差最大值误差最大值:0.0123 n 10 20 30 40 h 0.2 0.1 0.0667 0.05误差2 0.0123 0.0026 0.0011 5.9612e-004误差1 0.1059 0.0521 0.0342 0.025613/16欧拉方法欧拉方法,h=0.2时时误差最大值误差最大值:0.1059第13页，本讲稿共17页假设 yn=y(xn),称 Rn+1=y(xn+1)-yn+1为局部截断误差.即即由泰勒公式由泰勒公式Euler公式公式:yn+1=yn+hf(xn,yn)y(xn+1)yn+1=y(xn)yn+O(h2)=O(h2)14/16第14页，本讲稿共17页若局部截断误差为:O(h p+1)则称显式单步法具有则称显式单步法具有 p 阶精度阶精度 Euler方法的局部截断误差方法的局部截断误差:O(h 2)=O(h 1+1)所以所以,Euler公式是公式是 1 阶精度的方法阶精度的方法。验证验证 预测校正方法具有预测校正方法具有2阶精度阶精度15/16第15页，本讲稿共17页即 yn+1=yn+0.5hf(xn,yn)+f(xn+1,yn+hf(xn,yn)f(xn+1,yn+hf(xn,yn)=f(xn+h,yn+hf(xn,yn)=f(xn,yn)+hfxn+hf(xn,yn)fyn+O(h2)而而 y(xn+1)=y(xn)+hy(xn)+0.5h2y”(xn)+O(h3)0.5hf(xn,yn)+f(xn+1,yn+hf(xn,yn)=hy(xn)+0.5h2y”(xn)+O(h3)局部截断误差局部截断误差:y(xn+1)yn+1=y(xn)yn=O(h3)预测校正方法具有预测校正方法具有2阶精度阶精度16/16第16页，本讲稿共17页MATLAB解常微分方程初值问题命令求数值解命令求数值解命令:x,y=ode23(f,a,b,y0)f=inline(y-x.*y.2);x,y=ode23(f,0,2,1)解析解命令解析解命令:dsolve(eqn1,.)syms x ydsolve(Dy=y-x*y2,y(0)=1,x)ans=1/(x-1+2*exp(-x)解析解:第17页，本讲稿共17页
