一阶常微分方程的初值问题.ppt
一阶常微分方程的初值问题,Ordinary Differential Equations,ODE,一阶常微分方程的初值问题,一阶常微分方程的初值问题：节点：x1x2 xn步长 为常数,一阶常微分方程的初值问题,一 欧拉方法（折线法） 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),一阶常微分方程的初值问题,对于两个分量的一阶常微分方程组,用经典4阶 Runge-Kutta 法求解的格式为,一阶常微分方程的初值问题,n 级显式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, ,一阶常微分方程的初值问题,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, ,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 RKF45,一阶常微分方程的初值问题,Van 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);% define lags=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差分离散,一阶常微分方程的初值问题,一阶常微分方程的初值问题,一阶常微分方程的初值问题,一阶常微分方程的初值问题,Sod问题Sod问题是在激波管中充以两种介质，维持一定的压力差，用隔膜分开；当隔膜突然破裂后，形成间断面，研究其时间发展情况。Euler方程组：,
