一阶常微分方程的初值问题.pptx

Ordinary Differential EquationsODE一阶常微分方程的初值问题：节点：x1x2xn步长为常数一欧拉方法（折线法）yi+1=yi+h f(xi,yi)(i=0,1,n-1)优点：计算简单。缺点：一阶精度。二改进的欧拉方法改进的欧拉公式可改写为它每一步计算f(x,y)两次，截断误差为O(h3)精确解精确解：function t,y=Heun(ode,tspan,h,y0)t=(tspan(1):h:tspan(end);n=length(t);y=y0*ones(n,1);for i=2:n k1=feval(ode,t(i-1),y(i-1);k2=feval(ode,t(i),y(i-1)+h*k1);y(i)=y(i-1)+h*(k1+k2)/2;end三龙格库塔法(Runge-Kutta)欧拉公式可改写为它每一步计算f(xi,yi)一次，截断误差为O(h2)标准四阶龙格库塔公式每一步计算f(x,y)四次，截断误差为O(h5)01/21/21/201/210011/62/62/61/6对于两个分量的一阶常微分方程组对于两个分量的一阶常微分方程组用经典用经典4阶阶 Runge-Kutta 法求解的格式为法求解的格式为n 级显式级显式Runge-Kutta Runge-Kutta 方法的一般计算格式：方法的一般计算格式：其中Adams 外插公式（外插公式（Adams-Bashforth 公式）是一类公式）是一类 k+1 步步 k+1 阶显式方法阶显式方法三步法（三步法（k=2）,四步法（四步法（k=3）,Adams 内插公式（内插公式（Adams-Moulton 公式）是一类公式）是一类 k+1 步步 k+2 阶隐式方法阶隐式方法三步法（三步法（k=2）,Adams 预估预估-校正方法（校正方法（Adams-Bashforth-Moulton 公式）公式）一般取四步外插法与三步内插法结合。一般取四步外插法与三步内插法结合。#include#include#include#define TRUE 1main()int nstep_pr,j,k;float h,hh,k1,k2,k3,k4,t_old,t_limit,t_mid,t_new,t_pr,y,ya,yn;double fun();printf(n Fourth-Order Runge-Kutta Scheme n);while(TRUE)printf(Interval of t for printing?n);scanf(%f,&t_pr);printf(Number of steps in one printing interval?n);scanf(%d,&nstep_pr);printf(Maximum t?n);scanf(%f,&t_limit);y=1.0;/*Setting the initial value of the solution */h=t_pr/nstep_pr;printf(h=%g n,h);t_new=0;/*Time is initialized.*/hh=h/2;printf(-n);printf(t yn);printf(-n);printf(%12.5f%15.6e n,t_new,y);do for(j=1;j=nstep_pr;j+)t_old=t_new;t_new=t_new+h;yn=y;t_mid=t_old+hh;yn=y;k1=h*fun(yn,t_old);ya=yn+k1/2;k2=h*fun(ya,t_mid);ya=yn+k2/2;k3=h*fun(ya,t_mid);ya=yn+k3;k4=h*fun(ya,t_new);y=yn+(k1+k2*2+k3*2+k4)/6;printf(%12.5f%15.6e n,t_new,y);while(t_new=t_limit);printf(-n);printf(Maximum t limit exceeded n);printf(Type 1 to continue,or 0 to stop.n);scanf(%d,&k);if(k!=1)exit(0);double fun(y,t)float y,t;float fun_v;fun_v=-y;return(fun_v);四误差的控制我们常用事后估计法来估计误差，即从xi出发，用两种办法计算y(xi+1)的近似值。记为从xi出发以h为步长得到的y(xi+1)的近似值，记为从xi出发以h/2为步长跨两步得到的y(xi+1)的近似值。设给定精度为。如果不等式成立，则即为y(xi+1)的满足精度要求的近似值。自适应：自适应：使用使用2个不同的个不同的h。如果一个大的。如果一个大的h和一个和一个小的小的h得到的解相同，那么减小得到的解相同，那么减小h就没有意就没有意义了；相反如果两个解差别大，可以假设义了；相反如果两个解差别大，可以假设大大h值得到的解是不精确的。值得到的解是不精确的。使用相同的使用相同的h值，值，2种不同的算法。如果低种不同的算法。如果低精度算法与高精度算法的结果相同，则没精度算法与高精度算法的结果相同，则没有必要减小有必要减小h。Ode23 非刚性非刚性,单步法单步法,二三阶二三阶Runge-Kutta,精度低精度低Ode45非刚性非刚性,单步法单步法,四五阶四五阶Runge-Kutta,精度较高精度较高,最常用最常用Ode113非刚性非刚性,多步法多步法,采用可变阶采用可变阶(1-13)Adams PECE 算法算法,精度可高可低精度可高可低Ode15s 刚性刚性,多步法多步法,采用采用Gears(或或BDF)算法算法,精度精度中等中等.如果如果ode45很慢很慢,系统可能是刚性的系统可能是刚性的,可试此法可试此法Ode23s 刚性刚性,单步法单步法,采用采用2阶阶Rosenbrock法法,精度精度较低较低,可解决可解决ode15s 效果不好的刚性方程效果不好的刚性方程.Ode23t 适度刚性适度刚性,采用梯形法则采用梯形法则,适用于轻微刚性系统适用于轻微刚性系统,给出的解无数值衰减给出的解无数值衰减.Ode23tb 刚性刚性,TR-BDF2,即即R-K的第一级用梯形法则的第一级用梯形法则,第二级用第二级用Gear 法法.精度较低精度较低,对于误差允许范围比较差对于误差允许范围比较差的情况的情况,比比ode15s好好.Matlab 函数函数Matlabs ode23(Bogacki,Shampine)Runge-Kutta-Fehlberg方法方法(RKF45)4阶阶Runge-Kutta近似近似5阶阶Runge-Kutta近似近似局部误差估计局部误差估计Matlabs ode45 is a variation of RKF45Van der Pol:function xdot=vdpol(t,x)xdot(1)=x(2);xdot(2)=-(x(1)2-1)*x(2)-x(1);xdot=xdot;%VDPOL must return a column vector.%xdot=x(2);-(x(1)2-1)*x(2)-x(1);%xdot=0,1;-1,-(x(1)2-1)*x;t0=0;tf=20;x0=0;0.25;t,x=ode45(vdpol,t0,tf,x0);plot(t,x);figure(101)plot(x(:,1),x(:,2);Lorenz吸引子吸引子function xdot=lorenz(t,x)xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1*x;x0=0,0,eps;t,x=ode23(lorenz,0,100,x0);plot3(x(:,1),x(:,2),x(:,3);plot(x(:,1),x(:,2);function yp=examstiff(t,y)yp=-2,1;998,-998*y+2*sin(t);999*(cos(t)-sin(t);y0=2;3;tic,t,y=ode113(examstiff,0,10,y0);toctic,t,y=ode45(examstiff,0,10,y0);toctic,t,y=ode23(examstiff,0,10,y0);toctic,t,y=ode23s(examstiff,0,10,y0);toctic,t,y=ode15s(examstiff,0,10,y0);toctic,t,y=ode23t(examstiff,0,10,y0);toctic,t,y=ode23tb(examstiff,0,10,y0);toc刚性方程刚性方程向后差分方法（Gears method）隐式Runge-Kutta法function f=weissinger(t,y,yp)f=t*y2*yp3-y3*yp2+t*(t2+1)*yp-t2*y;t0=1;y0=sqrt(3/2);yp0=0;%guessy0,yp0=decic(weissinger,t0,y0,1,yp0,0);%求出自洽初值。保持y0不变t,y=ode15i(weissinger,1,10,y0,yp0);ytrue=sqrt(t.2+0.5);plot(t,ytrue,t,y,o)线性隐式线性隐式ODE:完全隐式完全隐式ODE(Matlab7):Weissinger方程方程:初值为时,解析解为 function yp=ddefun(t,y,Z)yp=zeros(2,1);%definelags=1,3yp(1)=Z(1,2)2+Z(2,1)2;yp(2)=y(1)+Z(2,1);function y=ddehist(t)y=zeros(2,1);y(1)=1;y(2)=t-2;lags=1,3;sol=dde23(ddefun,lags,ddehist,0,1);hold on;plot(sol.x,sol.y(1,:),b-);plot(sol.x,sol.y(2,:),r-);延迟微分方程延迟微分方程Sol=dde23(ddefun,lags,ddehist,tspan)初值:有限差分法有限差分法二阶线性边值问题二阶线性边值问题差分离散：差分离散：bvp4c线性边值问题的打靶法：二阶线性边值问题(11)的可以通过求解下面两个初值问题获得。原来边值问题的解可以表示为：非线性边值问题的打靶法(IVP1)(IVP2)符号计算符号计算y=dsolve(D2y=-a2*y,x)%求通解y=dsolve(D2y=-a2*y,y(0)=1,Dy(pi/a)=0,x)x,y=dsolve(Dx=4x-2y,Dy=2x-y,t)Pdetool求解区域,定义边界,网格划分,计算,画图,保存文件求解 边条 解析解 演示求解过程演示求解过程Stokes 问题问题Q1-P0有限元离散有限元离散Navier-Stokes 问题问题MAC差分离散差分离散物理问题的控制方程：前台阶流（前台阶流（A Mach 3 Wind Tunnel with a Step）模拟放置在风洞中的前台阶流动。风洞尺寸：宽1.0，长3.0，台阶高0.2，放置在距风洞左边界0.6个单位长度处。Sod问题问题Sod问题是在激波管中充以两种介质，维持一定的压力差，用隔膜分开；当隔膜突然破裂后，形成间断面，研究其时间发展情况。Euler方程组：A picture is worth a thousand words.-Anonymous Make it right before you make it faster.-Brain W.Kernighan,P.J.Plauger,The Elements of Programming Style(1978)