《MATLAB解微分方程》PPT课件.ppt

MATLABMATLABODEODE初值问题的数值解初值问题的数值解PDEPDE问题的数值解问题的数值解1精选课件ppt问题提出问题提出 倒葫芦形状容器壁上的刻度问题倒葫芦形状容器壁上的刻度问题.对于如图所示圆对于如图所示圆柱形状容器壁上的容积刻度柱形状容器壁上的容积刻度,可以利用圆柱体体积公式可以利用圆柱体体积公式其中直径其中直径D D为常数为常数.对于几何形状不是规则的容器对于几何形状不是规则的容器,比如倒葫芦形状容比如倒葫芦形状容器壁上如何标出刻度呢器壁上如何标出刻度呢?下表是经过测量得到部分下表是经过测量得到部分容器高度与直径的关系容器高度与直径的关系.2精选课件pptx1o根据上表的数据根据上表的数据,可以拟合出倒可以拟合出倒葫芦形状容器的图葫芦形状容器的图,建立如图所建立如图所示的坐标轴示的坐标轴,问题为如何根据任问题为如何根据任意高度意高度x x标出容器体积标出容器体积V V的刻度的刻度,由微元思想分析可知由微元思想分析可知H 0 0.2 0.4 0.6 0.8 1.0D 0.04 0.11 0.26 0.56 1.04 1.17其中其中x x表示高度表示高度,直径直径D D是高度是高度x x的函数的函数,记为记为D(x).D(x).x1o3精选课件ppt只要求解上述方程只要求解上述方程,就可求出体积就可求出体积V V与高度与高度x x之间的之间的函数关系函数关系,从而可标出容器壁上容积的刻度从而可标出容器壁上容积的刻度,但问题但问题是是函数函数D(x)D(x)无解析表达式无解析表达式,无法求出其解析解无法求出其解析解.因此因此,得到如下微分方程初值问题得到如下微分方程初值问题4精选课件ppt包含自变量、未知函数及未知函数的导数或微分的包含自变量、未知函数及未知函数的导数或微分的方程称为方程称为微分方程微分方程。在微分方程中在微分方程中,自变量的个数只有一个自变量的个数只有一个,称为称为常微常微分方程。分方程。自变量的个数为两个或两个以上的微分方程叫自变量的个数为两个或两个以上的微分方程叫偏微偏微分方程。分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的微分方程的阶数阶数。微分方程分类微分方程分类5精选课件ppt常微分方程：常微分方程：(1)(1)，(2)(2)式称为初值问题式称为初值问题.在实际应用中还经常需要求解常微分方程组：在实际应用中还经常需要求解常微分方程组：(3)(3)(3)(3)式称为边值问题。式称为边值问题。式称为边值问题。式称为边值问题。6精选课件ppt但能求但能求解析解解析解的常微分方程是有限的，大多数的常的常微分方程是有限的，大多数的常微分方程是给不出解析解的微分方程是给不出解析解的.这个一阶微分方程就不能用初等函数及其积分这个一阶微分方程就不能用初等函数及其积分来表达它的解。来表达它的解。例例 例例 的解的解的值仍需插值方法来计算的值仍需插值方法来计算.7精选课件ppt8精选课件ppt事实上，从事实上，从实际问题实际问题当中抽象出来的微分方程，通当中抽象出来的微分方程，通常主要依靠常主要依靠数值解法数值解法来解决。来解决。可以证明可以证明:如果函数在带形区域如果函数在带形区域 R=axb,R=axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数L(L(它与它与x,yx,y无关无关)使使 对对R R内任意两个内任意两个 都成立都成立,则方程的解则方程的解 在在 a,ba,b 上存在且唯一。上存在且唯一。在区间在区间a x ba x b上的数值解法上的数值解法。主要讨论一阶常微分方程初值问题主要讨论一阶常微分方程初值问题 9精选课件ppt微分方程数值方法的基本思想微分方程数值方法的基本思想 对常微分方程初值问题的数值解法，就是要算对常微分方程初值问题的数值解法，就是要算出精确解出精确解y(x)y(x)在区间在区间 a,ba,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值相邻两个节点的间距相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。假定长可以相等，也可以不等。假定h h为定数，称为定步长，为定数，称为定步长，这时节点可表示为这时节点可表示为10精选课件ppt数值解法需要把连续性的问题加以离散化，从数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。而求出离散节点的数值解。离散化离散化 对常微分方程数值解法的基本出发点就是对常微分方程数值解法的基本出发点就是离散化离散化。其数值解法的基本特点：采用其数值解法的基本特点：采用“步进式步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，即求解过程顺着节点排列的次序一步一步地向前推进，11精选课件ppt描述这类算法，要求给出用已知信息描述这类算法，要求给出用已知信息 计算计算 的递推公式。的递推公式。建立这类递推公式的基本方法是在这些节点上用建立这类递推公式的基本方法是在这些节点上用数数值积分、数值微分、泰勒展开值积分、数值微分、泰勒展开等离散化方法，对初等离散化方法，对初值问题值问题中的导数中的导数 进行不同的进行不同的离散化处理离散化处理。12精选课件ppt数值解和精确解数值解和精确解用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。用y(x)表示问题的准确解y(x0),y(x1),y(xN)表示解y(x)在节点x0,x1,xN处的准确值y0,y1,y N表示数值解，即问题的解y(x)在相应节点处的近似值。13精选课件ppt单步法和多步法单步法和多步法单步法：在计算yi+1 时只利用y i多步法：在计算yi+1 时不仅利用y i,还要利用 yi1,yi2,k步法：在计算yi+1 时要用到yi,yi1,yik+1显式格式可写成:yk+1=yk+hf(xk，yk;h）隐式格式:yk+1=yk+hf（xk，yk,yk+1;h）它每步求解yk+1需要解一个隐式方程。14精选课件ppt欧拉（欧拉（EulerEuler）方法）方法在x=x0 处，用差商代替导数：由由得得15精选课件ppt同理，在同理，在x=xx=xn n 处，用差商代替导数：处，用差商代替导数：由由得得若记若记则上式可记为则上式可记为此此即即为为求求解解初初值值问问题题的的EulerEuler方方法法，又又称称显显式式EulerEuler方法。方法。16精选课件ppt例例:用用EulerEuler方法求解常微分方程初值问题方法求解常微分方程初值问题并将数值解和该问题的解析解比较。并将数值解和该问题的解析解比较。解：解：EulerEuler方法的具体格式：方法的具体格式：17精选课件ppth=0.2;y(1)=0.2;x=0.2:h:3;h=0.2;y(1)=0.2;x=0.2:h:3;for n=1:14for n=1:14 xn=x(n);yn=y(n);xn=x(n);yn=y(n);y(n+1)=yn+h*(yn/xn-2*yn*yn);y(n+1)=yn+h*(yn/xn-2*yn*yn);endendx0=0.2:h:3;y0=x0./(1+x0.2);x0=0.2:h:3;y0=x0./(1+x0.2);plot(x0,y0,x,y,x,y,o)plot(x0,y0,x,y,x,y,o)程序实现程序实现18精选课件pptxn y(xn)yn yn-y(xn)0.00000.20.19230.20000.00770.40.34480.38400.03920.60.44120.51700.07580.80.48780.58240.09461.00.50000.59240.09241.20.49180.57050.07871.40.47300.53540.0624h h=0.2,=0.2,x xn n=nhnh,(,(n n=0,1,2,15),=0,1,2,15),f f(x x,y y)=)=y y/x x 2 2y y2 2 计算中取计算中取f f(0,0)=1.(0,0)=1.计算结果如下：计算结果如下：19精选课件pptxn y(xn)yn yn-y(xn)1.60.44940.49720.04781.80.42450.46050.03592.00.40000.42680.02682.20.37670.39660.01992.40.35500.36980.01472.60.33510.34590.01082.80.31670.32460.00793.00.30000.30570.0057由由表表中中数数据据可可以以看看到到，微微分分方方程程初初值值问问题题的的数数值值解解和和解解析析解解的的误误差差一一般般在在小小数数点点后后第第二二位位或或第第三三位位小小数数上上，说说明明EulerEuler方法的精度是比较差的。方法的精度是比较差的。20精选课件ppt二阶二阶Runge-KuttaRunge-Kutta方法方法其中其中其中其中 c c1 1,c c2 2,2 2,21 21 待定。待定。待定。待定。上式的上式的局部截断误差局部截断误差为：为：21精选课件ppt即即c c1 1=1-1-a a，2 2=21 21 =1/(21/(2a a)方程组解不唯一，可令方程组解不唯一，可令c c2 2=a=a 0 0 ，则，则 满足上述条件的公式都为满足上述条件的公式都为2 2阶阶R-KR-K公式公式。22精选课件ppt例例 蛇形曲线蛇形曲线的初值问题的初值问题令令f(x,y)=y/x 2yf(x,y)=y/x 2y2 2,取取 f(0,0)=1,f(0,0)=1,h=0.2,xh=0.2,xn n=hn ,(n=1,2,15)=hn ,(n=1,2,15)2 2阶龙格阶龙格-库塔公式计算格式：库塔公式计算格式：k k1 1=y=yn n/x/xn n 2y 2yn n 2 2,k k2 2=(y=(yn n+hk+hk1 1)/(x)/(xn n+h)2+h)2(y(yn n+hk+hk1 1)2 2 y yn+1n+1=y=yn n+0.5h k+0.5h k1 1+k+k2 2 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);for n=1:14 k1=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-KuttaRunge-Kutta方法方法25精选课件ppt function ydot=harmonic(t,y)ydot=y(2);-y(1)y=inline(0 1;-1 0*y,t,y);System of Equations26精选课件ppt function ydot=twobody(t,y)r=sqrt(y(1)2+y(2)2);ydot=y(3);y(4);-y(1)/r3;-y(2)/r3;Two Body Problem27精选课件pptLinearized Differential Equations28精选课件pptJ的特征值是解增长解衰减解振荡29精选课件ppt 基于龙格库塔法，基于龙格库塔法，MATLABMATLAB求常微分方程数值解求常微分方程数值解的函数的函数，一般调用格式为：一般调用格式为：t,y=ode23(fname,tspan,y0)t,y=ode23(fname,tspan,y0)t,y=ode45(fname,tspan,y0)t,y=ode45(fname,tspan,y0)其中其中fnamefname是定义是定义f(t,y)f(t,y)的函数文件名的函数文件名，该函数文件，该函数文件必须返回一个列向量。必须返回一个列向量。tspantspan形式为形式为t0,tf,t0,tf,表示表示求解区间。求解区间。y0y0是初始状态列向量。是初始状态列向量。t t和和y y分别给出分别给出时间向量和相应的状态向量。时间向量和相应的状态向量。MATLABMATLAB求常微分方程数值解的函数求常微分方程数值解的函数30精选课件ppt ode23:Bogacki,Shampine(1989)和Shampine(1994),”23”表示用两算法：一个2阶，一个3阶Bogacki,P.and LF Shampine,A 3(2)pair of Runge-Kutta formulas,Appl.Math.Letters,Vol.2,1989,pp 1-9.BS23 algorithm31精选课件ppt F=inline(y(2);-y(1),t,y)ode23(F,0 2*pi,1;0)opts=odeset(reltol,1.e-4,abstol,1.e-6,outputfcn)Examples32精选课件ppt ode23(twobody,0 2*pi,1;0;0;1);01234567-1.5-1-0.500.511.5Examples33精选课件ppt y0=1;0;0;3;ode23(twobody,0 2*pi,y0);01234567-202468101214161834精选课件ppt y0=1;0;0;3;t,y=ode23(twobody,0 2*pi,y0);plot(y(:,1),y(:,2);axis equal35精选课件ppt y0=1;0;0;3;t,y=ode23(twobody,0 2*pi,y0);plot(y(:,1)plot(y(:,2)36精选课件ppt A problem is stiff if the solution being sought is varying slowly,but there are nearby solutions that very rapidly,so the numerical method must take small steps to obtain satisfactory results.A model of flame propagation（火焰燃烧）:y是球的半径，y2和y3与球的表面积和体积有关想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小（由于进入球内氧气和消耗的氧气平衡）Stiff Problem(刚性问题)37精选课件ppt010203040506070809010000.20.40.60.811.21.4eta=0.02;sym y;F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);38精选课件ppteta=0.00002;sym y;F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);39精选课件ppt012345678910 x 10400.20.40.60.811.21.4ode23seta=0.00002;ode23s(inline(y2-y3,t,y),0 2/eta,eta);40精选课件ppt例例 蛇形曲线的常微分方程初值问题蛇形曲线的常微分方程初值问题 MATLABMATLAB数值求解命令数值求解命令F=inline(1./(1+x.2)-2*y.2)F=inline(1./(1+x.2)-2*y.2)；ode23(F,0,6,0)ode23(F,0,6,0)输出结果为图形输出结果为图形 41精选课件pptT,y=ode23(f,0,6,0)T,y=ode23(f,0,6,0)将将得到自变量和函数的离散数得到自变量和函数的离散数据据 T=0 0.0001 0.0005 0.0025 0.0125 0.0496 0.1085 0.1863 0.2837 0.4091 0.5991 0.8513 1.0567 1.2680 1.5110 1.8050 2.1788 2.6842 3.2842 3.8842 4.4842 5.0842 5.6842 6.0000y=0 0.0001 0.0005 0.0025 0.0125 0.0495 0.1073 0.1800 0.2626 0.3505 0.4411 0.4944 0.5000 0.4868 0.4607 0.4242 0.3793 0.3270 0.2783 0.2411 0.2122 0.1891 0.1705 0.162042精选课件ppt例例 洛伦兹模型洛伦兹模型由如下常微分方程组描述由如下常微分方程组描述 取取 =8/38/3，=1010，=2828。初值初值：x x(0)=0(0)=0，y y(0)=0(0)=0，z z(0)=0.01(0)=0.01。利用利用MATLABMATLAB求解常微分方程数值解命令计算出求解常微分方程数值解命令计算出t t00，8080内内，三个未知函数的数据值，并绘出相空间在三个未知函数的数据值，并绘出相空间在Y-XY-X平面的投影曲线平面的投影曲线 43精选课件ppt气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应蝴蝶效应。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。44精选课件ppt天气预报的准确性：http:/.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLorenz现象的数学：http:/www.sciscape.org/news_detail.php?news_id=225分形艺术电子版：http:/ y y1 1，y y2 2，y y3 3=x x，y y，z z，创建创建MATLABMATLAB函函数文件如下数文件如下 function z=flo(t,y)function z=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3);z(2,:)=-10*(y(2)-y(3);z(3,:)=-y(1).*y(2)+28*y(2)-y(3);z(3,:)=-y(1).*y(2)+28*y(2)-y(3);用用MATLABMATLAB命令求解并绘出命令求解并绘出Y-XY-X平面的投影图平面的投影图 y0=0;0;0.01;y0=0;0;0.01;x,y=ode23(flo,0,80,y0)x,y=ode23(flo,0,80,y0)plot(y(:,2),y(:,1)plot(y(:,2),y(:,1)46精选课件pptplot(y(:,3),y(:,1)plot(y(:,3),y(:,1)47精选课件pptplot3(y(:,1),y(:,2),y(:,3)48精选课件ppt y0=30;0;-40;plot(y(:,i)y0=30;0;-40;plot(y(:,i)49精选课件ppt50精选课件ppt51精选课件ppt非刚性系统：ode45(Runge-Kutta45)ode23(Runge-Kuatta23)ode113(Adams-Bashforth-Moulton PECE)多步方法刚性系统：ode15s(Gear方法),多步方法 ode23s(二阶modified Rosenbroack formula)，单步 ode23t(trapezoidal rule),solve DAEs ode23tb(TR-BDF2)low order methodMatlabs ODE Solvers52精选课件pptLaplacian 算子：Poisson方程（elliptic）：Laplacian 算子的特征值问题:Heat equation(parabolic):Wave equation(hyperbolic):PDE Model53精选课件ppt五点离散 Poisson方程离散：特征值问题:Finete Difference Methods54精选课件ppt热方程：波动方程：55精选课件ppt-4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1-4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0-4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1-4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0-4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1-4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0-4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1-4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1-4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1-4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1-4 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0-4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1-4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1-4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1-4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1-4Matrixs56精选课件ppt椭圆方程：特征值方程：热方程：波动方程：离散方程57精选课件ppt 例 用古典显式格式求解抛物型方程初始条件为边界条件为取步长x=h=0.2,t=k=0.02。58精选课件ppt 解 r=k/h2=0.02/0.22=0.5,古典显式格式为 0,0 0.2,0 0.4,0 1,0 1,00,0.20,0.020,0.021,0.21,0.0259精选课件pptfunction u=gu_dian(f,a,b,c1,c2,m,n)%输入初值和Uh=a/(m-1);k=b/(n-1);r=k/h2;U=zeros(n,m);%赋边界条件U(2:n,1)=c1;U(2:n,m)=c2;function y=fg(x)y=4.*(x-x.2);60精选课件ppt%赋初始条件U(1,1:m)=fg(0:h:h*(m-1);%计算内点上u的数值解Ufor i=2:n for j=2:(m-1)U(i,j)=(1-2*r)*U(i-1,j)+r*(U(i-1,j-1)+U(i-1,j+1);endend61精选课件ppt%gu_dianl1.m 步长h=0.20,k=0.02,r=k/h2=0.5a=1;b=0.20;c1=0;c2=0;m=6;n=11;U=gu_dian(fg,a,b,c1,c2,m,n)x=0:0.2:a;y=0:0.02:b;X,Y=meshgrid(x,y);surf(X,Y,U)%输入U后再画图62精选课件ppt63精选课件ppt64精选课件pptPDETOOLPDETOOL有限元方法有限元方法65精选课件ppt