数值分析常微分方程初值问题幻灯片.ppt

《数值分析常微分方程初值问题幻灯片.ppt》由会员分享，可在线阅读，更多相关《数值分析常微分方程初值问题幻灯片.ppt（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数值分析常微分方程初值问题第1页，共17页，编辑于2022年，星期六例例1.罗伦茨模型罗伦茨模型求解区域求解区域:0 t 80初始条件：初始条件：x(0)=0,y(0)=0,z(0)=0.0001取取 =8/3,=28,=102/16第2页，共17页，编辑于2022年，星期六例例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页，编辑于2022年，星期六例例3.人口模型人口模型(Logistic):P=r

2、P(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页，编辑于2022年，星期六常微分方程初值问题常微分方程初值问题:例例 数值方法：数值方法：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页，编辑于2022年，星期六连续问题连续问题:离散问题离散问题:n=1,2,3,6/16第6页，共17

3、页，编辑于2022年，星期六求解常微分方程初值问题的求解常微分方程初值问题的Euler方法方法 7/16第7页，共17页，编辑于2022年，星期六例例 取步长取步长 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页，编辑于2022年，星期六n=inp

4、ut(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页，编辑于2022年，星期六y=f(x,y)梯形公式梯形公式:左矩形公式左矩形公式用数值积分方

5、法离散化常微分方程用数值积分方法离散化常微分方程10/16第10页，共17页，编辑于2022年，星期六算法算法:k1=f(xn,yn),k2=f(xn+1,yn+h k1)由梯形公式推出的预由梯形公式推出的预-校方法校方法:11/16第11页，共17页，编辑于2022年，星期六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*

6、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页，编辑于2022年，星期六预预-校方法校方法,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-

7、004误差误差1 0.1059 0.0521 0.0342 0.025613/16欧拉方法欧拉方法,h=0.2时时误差最大值误差最大值:0.1059第13页，共17页，编辑于2022年，星期六假设假设 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页，编辑于2022年，星期六若局部截断误差为若局部截断误差为:O(h p+1)则称显式单步法具有则称显式单步法具有 p 阶精度阶精度 Euler

8、方法的局部截断误差方法的局部截断误差:O(h 2)=O(h 1+1)所以所以,Euler公式是公式是 1 阶精度的方法阶精度的方法。验证验证 预测校正方法具有预测校正方法具有2阶精度阶精度15/16第15页，共17页，编辑于2022年，星期六即即 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+

9、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页，编辑于2022年，星期六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页，编辑于2022年，星期六