Ordinary Differential EquationsODE 一阶常微分方程的初值问题： 节点：x1 #include #include #define TRUE 1 main() 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“, printf( “Number of steps in one printing interval?n“ );scanf( “%d“, printf( “Maximum t?n“ );scanf( “%f“, 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 );dofor( 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“, 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);toc tic,t,y = ode45(examstiff,0,10,y0);toc tic,t,y = ode23(examstiff,0,10,y0);toc tic,t,y = ode23s(examstiff,0,10,y0);toc tic,t,y = ode15s(examstiff,0,10,y0);toc tic,t,y = ode23t(examstiff,0,10,y0);toc tic,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; % guess y0,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,3 yp(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)