一阶常微分方程的初值问题课件

OrdinaryDifferentialEquationsODE一阶常微分方程的初值问题：节点：x1<x2<…<xn步长为常数一欧拉方法（折线法）

yi+1=yi+hf(xi,yi)(i=0,1,…,n-1)优点：计算简单。缺点：一阶精度。二改进的欧拉方法精确解：function[t,y]=Heun(ode,tspan,h,y0)t=(tspan(1):h:tspan(end))';n=length(t);y=y0*ones(n,1);fori=2:nk1=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/6n

级显式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<stdio.h>#include<stdlib.h>#include<math.h>#defineTRUE1main(){intnstep_pr,j,k;floath,hh,k1,k2,k3,k4,t_old,t_limit,t_mid,t_new,t_pr,y,ya,yn;doublefun();printf("\nFourth-OrderRunge-KuttaScheme\n");while(TRUE){printf("Intervaloftforprinting?\n"); scanf("%f",&t_pr);printf("Numberofstepsinoneprintinginterval?\n"); scanf("%d",&nstep_pr);printf("Maximumt?\n"); scanf("%f",&t_limit);y=1.0;/*Settingtheinitialvalueofthesolution*/h=t_pr/nstep_pr;printf("h=%g\n",h);t_new=0;/*Timeisinitialized.*/hh=h/2;printf("--------------------------------------\n");printf("ty\n");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("Maximumtlimitexceeded\n");printf("Type1tocontinue,or0tostop.\n");scanf("%d",&k);if(k!=1)exit(0);}}doublefun(y,t)floaty,t;{floatfun_v;fun_v=-y;return(fun_v);}自适应：使用2个不同的h。如果一个大的h和一个小的h得到的解相同，那么减小h就没有意义了；相反如果两个解差别大，可以假设大h值得到的解是不精确的。使用相同的h值，2种不同的算法。如果低精度算法与高精度算法的结果相同，则没有必要减小h。Ode23非刚性,单步法,二三阶Runge-Kutta,精度低Ode45非刚性,单步法,四五阶Runge-Kutta,精度较高,最常用Ode113非刚性,多步法,采用可变阶(1-13)AdamsPECE算法,精度可高可低Ode15s刚性,多步法,采用Gear’s(或BDF)算法,精度中等.如果ode45很慢,系统可能是刚性的,可试此法Ode23s刚性,单步法,采用2阶Rosenbrock法,精度较低,可解决ode15s效果不好的刚性方程.Ode23t适度刚性,采用梯形法则,适用于轻微刚性系统,给出的解无数值衰减.Ode23tb刚性,TR-BDF2,即R-K的第一级用梯形法则,第二级用Gear法.精度较低,对于误差允许范围比较差的情况,比ode15s好.Matlab函数Matlab’sode23(Bogacki,Shampine)VanderPol:functionxdot=vdpol(t,x)xdot(1)=x(2);xdot(2)=-(x(1)^2-1)*x(2)-x(1);xdot=xdot';%VDPOLmustreturnacolumnvector. %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吸引子functionxdot=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));functionf=weissinger(t,y,yp)f=t*y^2*yp^3-y^3*yp^2+t*(t^2+1)*yp-t^2*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方程:初值为时,解析解为functionyp=ddefun(t,y,Z)yp=zeros(2,1);%definelags=[1,3]yp(1)=Z(1,2)^2+Z(2,1)^2;yp(2)=y(1)+Z(2,1);functiony=ddehist(t)y=zeros(2,1);y(1)=1;y(2)=t-2;lags=[1,3];sol=dde23(@ddefun,lags,@ddehist,[0,1]);holdon;plot(sol.x,sol.y(1,:),'b-');plot(sol.x,sol.y(2,:),'r-');延迟微分方程Sol=dde23(ddefun,lags,ddehist,tspan)初值

:有限差分法二阶线性边值问题差分离散：bvp4c符号计算y=dsolve('D2y=-a^2*y','x')%求通解y=dsolve('D2y=-a^2*y','y(0)=1','Dy(pi/a)=0','x')[x,y]=dsolve('Dx=4x-2y','Dy=2x-y','t')Pdetool求解区域,定义边界,网格划分,计算,画图,保存文件求解边条解析解演示求解过程Navier-Stokes问题MAC差分离散物理问题的控制方程：前台阶流（AMach3WindTunnelwithaStep）模拟放置在风洞中的前台阶流动。风洞尺寸：宽1.0，长3.0，台阶高0.2，放置在距风洞左边界0.6个单位长度处。Sod问题Sod问题是在激波管中充以两种介质，维持一定的压力差，用隔膜分开；当隔膜突然破裂后，形成间断面，研究其时间发展情况。Euler方程组：Apictureiswortha