《一阶常微方程的初值问题ppt课件》由会员分享，可在线阅读，更多相关《一阶常微方程的初值问题ppt课件（37页珍藏版）》请在金锄头文库上搜索。
1、Ordinary Differential EquationsODE一阶常微分方程的初值问题：节点：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
2、-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 方法的普通计算格式：方法的
3、普通计算格式：其中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
4、, 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( M
5、aximum 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_o
6、ld = 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
7、= 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)的近似值，记 为从x
8、i出发以 h/2 为步长跨两步得到的y(xi+1)的近似值。设给定精度为。假设不等式 成立，那么 即为y(xi+1)的满足精度要求的近似值。自顺应：自顺应：运用运用2个不同的个不同的h。假设一个。假设一个大的大的h和一个小的和一个小的h得到的得到的解一样，那么减小解一样，那么减小h就没有就没有意义了；相反假设两个解差意义了；相反假设两个解差别大，可以假设大别大，可以假设大h值得到值得到的解是不准确的。的解是不准确的。运用一样的运用一样的h值，值，2种不同的种不同的算法。假设低精度算法与高算法。假设低精度算法与高精度算法的结果一样，那么精度算法的结果一样，那么没有必要减小没有必要减小h。Ode2
9、3 非非刚性性, 单步法步法, 二三二三阶Runge-Kutta,精度低精度低Ode45非非刚性性, 单步法步法, 四五四五阶Runge-Kutta,精度精度较高高,最常用最常用Ode113非非刚性性, 多步法多步法, 采用可采用可变阶(1-13)Adams PECE 算法算法, 精度可高可低精度可高可低Ode15s 刚性性, 多步法多步法,采用采用Gears (或或BDF)算法算法, 精度精度中等中等. 假假设ode45很慢很慢, 系系统能能够是是刚性的性的,可可试此法此法Ode23s 刚性性, 单步法步法, 采用采用2阶Rosenbrock法法, 精度精度较低低, 可可处理理ode15s
10、效果不好的效果不好的刚性方程性方程.Ode23t 适度适度刚性性, 采用梯形法那么采用梯形法那么,适用于适用于细微微刚性系性系统,给出的解无数出的解无数值衰减衰减.Ode23tb 刚性性, TR-BDF2, 即即R-K的第一的第一级用梯形法那用梯形法那么么,第二第二级用用Gear 法法. 精度精度较低低, 对于于误差允差允许范范围比比较差的情况差的情况,比比ode15s好好.Matlab 函数函数Matlabs ode23(Bogacki, Shampine)Runge-Kutta-Fehlberg方法方法(RKF45)4阶阶Runge-Kutta近似近似5阶阶Runge-Kutta近似近似部
11、分误差估计部分误差估计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(
12、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*(
13、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刚性方程刚
14、性方程向后差分方法Gearsmethod隐式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
15、: 完全完全隐式式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(
16、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,
17、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)
