《MATLAB解微分方程》课件

MATLABODE初值问题的数值解PDE问题的数值解1精选课件ppt问题提出倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数.对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下表是经过测量得到部分容器高度与直径的关系.2精选课件pptx1o根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴,问题为如何根据任意高度x标出容器体积V的刻度,由微元思想分析可知H00.20.40.60.81.0D0.040.110.260.561.041.17其中x表示高度,直径D是高度x的函数,记为D(x).x1o3精选课件ppt只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是函数D(x)无解析表达式,无法求出其解析解.

因此,得到如下微分方程初值问题4精选课件ppt包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。微分方程分类5精选课件ppt常微分方程：(1)，(2)式称为初值问题.在实际应用中还经常需要求解常微分方程组：(3)式称为边值问题。6精选课件ppt但能求解析解的常微分方程是有限的，大多数的常微分方程是给不出解析解的.这个一阶微分方程就不能用初等函数及其积分来表达它的解。

例例

的解的值仍需插值方法来计算.7精选课件ppt8精选课件ppt事实上，从实际问题当中抽象出来的微分方程，通常主要依靠数值解法来解决。可以证明:如果函数在带形区域R=a≤x≤b,-∞＜y＜∞}内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使

对R内任意两个都成立,则方程的解在

a,b

上存在且唯一。

在区间a≤x≤b上的数值解法。

主要讨论一阶常微分方程初值问题9精选课件ppt微分方程数值方法的基本思想对常微分方程初值问题的数值解法，就是要算出精确解y(x)在区间

a,b

上的一系列离散节点

处的函数值相邻两个节点的间距称为步长，步长可以相等，也可以不等。假定h为定数，称为定步长，这时节点可表示为10精选课件ppt数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。离散化对常微分方程数值解法的基本出发点就是离散化。其数值解法的基本特点：采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，11精选课件ppt描述这类算法，要求给出用已知信息计算的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数进行不同的离散化处理。12精选课件ppt数值解和精确解用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。用y(x)表示问题的准确解y(x0),y(x1),y(xN)表示解y(x)在节点x0,x1,…,xN处的准确值y0,y1,…,yN表示数值解，即问题的解y(x)在相应节点处的近似值。13精选课件ppt单步法和多步法单步法：在计算yi+1时只利用yi多步法：在计算yi+1时不仅利用yi,还要利用yi−1,yi−2,…,k步法：在计算yi+1时要用到yi,yi−1,…,yi−k+1显式格式可写成:yk+1=yk+hΦf(xk，yk;h）隐式格式:yk+1=yk+hΦf（xk，yk,yk+1;h）它每步求解yk+1需要解一个隐式方程。14精选课件ppt欧拉（Euler）方法在x=x0处，用差商代替导数：由得15精选课件ppt同理，在x=xn

处，用差商代替导数：由得若记则上式可记为此即为求解初值问题的Euler方法，又称显式Euler方法。16精选课件ppt例:用Euler方法求解常微分方程初值问题并将数值解和该问题的解析解比较。解：Euler方法的具体格式：17精选课件ppth=0.2;y(1)=0.2;x=0.2:h:3;forn=1:14xn=x(n);yn=y(n);y(n+1)=yn+h*(yn/xn-2*yn*yn);endx0=0.2:h:3;y0=x0./(1+x0.^2);plot(x0,y0,x,y,x,y,'o')程序实现18精选课件ppt

xn

y(xn) yn

yn-y(xn) 0.0 0 0 00.2 0.1923 0.2000 0.00770.4 0.3448 0.3840 0.03920.6 0.4412 0.5170 0.07580.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.09241.2 0.4918 0.5705 0.07871.4 0.4730 0.5354 0.0624h=0.2,xn=nh,(n=0,1,2…,15),f(x,y)=y/x–2y2计算中取f(0,0)=1.计算结果如下：19精选课件pptxn

y(xn) yn

yn-y(xn)1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.03592.0 0.4000 0.4268 0.02682.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.01472.6 0.3351 0.3459 0.01082.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小数上，说明Euler方法的精度是比较差的。20精选课件ppt二阶Runge-Kutta方法其中c1,c2,

2,

21待定。上式的局部截断误差为：21精选课件ppt即c1=1-a，

2=

21=1/(2a)方程组解不唯一，可令c2=a

0

，则

满足上述条件的公式都为2阶R-K公式。22精选课件ppt例

蛇形曲线的初值问题令f(x,y)=y/x–2y2,取f(0,0)=1,h=0.2,xn=hn,(n=1,2,…,15)2阶龙格-库塔公式计算格式：k1=yn/xn–2yn

2,k2=(yn+hk1)/(xn+h)–2(yn+hk1)2yn+1=yn+0.5h[k1+k2]23精选课件pptx0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2;y(1)=y0+.5*h*(k1+k2);forn=1:14k1=y(n)/x(n)-2*y(n)^2;k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2;y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.^2);plot(x,y,'o',x,y1)24精选课件ppt常用的一个公式为四阶Runge-Kutta方法25精选课件pptfunctionydot=harmonic(t,y)ydot=[y(2);-y(1)]y=inline(‘[01;-10]*y’,’t’,’y’);SystemofEquations26精选课件pptfunctionydot=twobody(t,y)r=sqrt(y(1)^2+y(2)^2);ydot=[y(3);y(4);-y(1)/r^3;-y(2)/r^3];

TwoBodyProblem27精选课件pptLinearizedDifferentialEquations28精选课件pptJ的特征值是解增长解衰减解振荡29精选课件ppt基于龙格－库塔法，MATLAB求常微分方程数值解的函数，一般调用格式为：[t,y]=ode23('fname',tspan,y0)[t,y]=ode45('fname',tspan,y0)其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。MATLAB求常微分方程数值解的函数30精选课件pptode23:Bogacki,Shampine(1989)和Shampine(1994),”23”表示用两算法：一个2阶，一个3阶Bogacki,P.andLFShampine,"A3(2)pairofRunge-Kuttaformulas,"Appl.

Math.Letters,Vol.2,1989,pp1-9.

BS23algorithm31精选课件pptF=inline('[y(2);-y(1)]','t','y')ode23(F,[02*pi],[1;0])opts=odeset(‘reltol’,1.e-4,’abstol’,1.e-6,’outputfcn’)Examples32精选课件pptode23(@twobody,[02*pi],[1;0;0;1]);01234567-1.5-1-0.500.511.5Examples33精选课件ppty0=[1;0;0;3];ode23(@twobody,[02*pi],y0);01234567-202468101214161834精选课件ppty0=[1;0;0;3];[t,y]=ode23(@twobody,[02*pi],y0);plot(y(:,1),y(:,2));axisequal35精选课件ppty0=[1;0;0;3];[t,y]=ode23(@twobody,[02*pi],y0);plot(y(:,1))plot(y(:,2))36精选课件pptAproblemisstiffifthesolutionbeingsoughtisvaryingslowly,buttherearenearbysolutionsthatveryrapidly,sothenumericalmethodmusttakesmallstepstoobtainsatisfactoryresults.Amodelofflamepropagation（火焰燃烧）:y是球的半径，y^2和y^3与球的表面积和体积有关想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小（由于进入球内氧气和消耗的氧气平衡）StiffProblem(刚性问题)37精选课件ppt010203040506070809010000.20.40.60.811.21.4eta=0.02;symy;F=inline(‘y^2-y^3’,’t’,’y’);ode23(F,[02/eta],eta);38精选课件ppteta=0.00002;symy;F=inline(‘y^2-y^3’,’t’,’y’);ode23(F,[02/eta],eta);39精选课件ppt012345678910x10400.20.40.60.811.21.4ode23seta=0.00002;ode23s(inline('y^2-y^3','t','y'),[02/eta],eta);40精选课件ppt例蛇形曲线的常微分方程初值问题MATLAB数值求解命令F=inline('1./(1+x.^2)-2*y.^2')；ode23(F,[0,6],0)输出结果为图形41精选课件ppt[T,y]=ode23(f,[0,6],0)将得到自变量和函数的离散数据T=00.00010.00050.00250.01250.04960.10850.18630.28370.40910.59910.85131.05671.26801.51101.80502.17882.68423.28423.88424.48425.08425.68426.0000y=00.00010.00050.00250.01250.04950.10730.18000.26260.35050.44110.49440.50000.48680.46070.42420.37930.32700.27830.24110.21220.18910.17050.162042精选课件ppt例洛伦兹模型由如下常微分方程组描述取

=8/3，

=10，

=28。初值：x(0)=0，y(0)=0，z(0)=0.01。利用MATLAB求解常微分方程数值解命令计算出t∈[0，80]内，三个未知函数的数据值，并绘出相空间在Y-X平面的投影曲线43精选课件ppt气象学家Lorenz提出一篇论文，名叫「一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？」论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做「蝴蝶效应」。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。44精选课件ppt天气预报的准确性：.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLorenz现象的数学：/news_detail.php?news_id=225分形艺术电子版：/personal/huajie/fractalart/html/notice.htm混沌理论：/~yzhao/doc/chaos.html45精选课件ppt记向量[y1，y2，y3]=[x，y，z]，创建MATLAB函数文件如下functionz=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3));z(3,:)=-y(1).*y(2)+28*y(2)-y(3);用MATLAB命令求解并绘出Y-X平面的投影图

y0=[0;0;0.01];[x,y]=ode23(@flo,[0,80],y0)plot(y(:,2),y(:,1))46精选课件pptplot(y(:,3),y(:,1))47精选课件pptplot3(y(:,1),y(:,2),y(:,3))48精选课件ppt

y0=[30;0;-40];plot(y(:,i))49精选课件ppt50精选课件ppt51精选课件ppt非刚性系统：ode45(Runge-Kutta45)ode23(Runge-Kuatta23)ode113(Adams-Bashforth-MoultonPECE)多步方法刚性系统：ode15s(Gear方法),多步方法ode23s(二阶modifiedRosenbroackformula)，单步ode23t(trapezoidalrule),solveDAEsode23tb(TR-BDF2)lowordermethodMatlab'sODESolvers52精选课件pptLaplacian算子：Poisson方程（elliptic）：Laplacian算子的特征值问题:Heatequation(parabolic):Waveequation(hyperbolic):PDEModel53精选课件ppt五点离散Poisson方程离散：特征值问题:FineteDifferenceMethods54精选课件ppt热方程：波动方程：55精选课件ppt-41100000000000001-40100000000000010-41100000000000011-40100000000000010-41100000000000011-40100000000000010-41000100000000011-41000100000000001-41000100000000001-41000100000000001-40000100000000000-41000000000000001-41000000000000001-4100000