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

《《解微分方程》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《解微分方程》PPT课件.ppt（65页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、MATLABMATLABODEODE初值问题的数值解初值问题的数值解PDEPDE问题的数值解问题的数值解1问题提出问题提出 倒葫芦形状容器壁上的刻度问题倒葫芦形状容器壁上的刻度问题.对于如图所示圆对于如图所示圆柱形状容器壁上的容积刻度柱形状容器壁上的容积刻度,可以利用圆柱体体积公式可以利用圆柱体体积公式其中直径其中直径D D为常数为常数.对于几何形状不是规则的容器对于几何形状不是规则的容器,比如倒葫芦形状容比如倒葫芦形状容器壁上如何标出刻度呢器壁上如何标出刻度呢?下表是经过测量得到部分下表是经过测量得到部分容器高度与直径的关系容器高度与直径的关系.2x1o根据上表的数据根据上表的数据,可以拟合

2、出倒可以拟合出倒葫芦形状容器的图葫芦形状容器的图,建立如图所建立如图所示的坐标轴示的坐标轴,问题为如何根据任问题为如何根据任意高度意高度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只要求解上述方程只要求解上述方程,就可求出体积就可求出体积V V与高度与高度x x之间的之间的函数关系函数关系,从而可标出容器壁上容积的刻度从而可标出容器壁

3、上容积的刻度,但问题但问题是是函数函数D(x)D(x)无解析表达式无解析表达式,无法求出其解析解无法求出其解析解.因此因此,得到如下微分方程初值问题得到如下微分方程初值问题4包含自变量、未知函数及未知函数的导数或微分的包含自变量、未知函数及未知函数的导数或微分的方程称为方程称为微分方程微分方程。在微分方程中在微分方程中,自变量的个数只有一个自变量的个数只有一个,称为称为常微常微分方程。分方程。自变量的个数为两个或两个以上的微分方程叫自变量的个数为两个或两个以上的微分方程叫偏微偏微分方程。分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的微

4、分方程的阶数阶数。微分方程分类微分方程分类5常微分方程常微分方程：(1)(1)，(2)(2)式称为初值问题式称为初值问题.在实际应用中还经常需要求解常微分方程组：在实际应用中还经常需要求解常微分方程组：(3)(3)(3)(3)式称为边值问题。式称为边值问题。式称为边值问题。式称为边值问题。6但能求但能求解析解解析解的常微分方程是有限的，大多数的常的常微分方程是有限的，大多数的常微分方程是给不出解析解的微分方程是给不出解析解的.这个一阶微分方程就不能用初等函数及其积这个一阶微分方程就不能用初等函数及其积分来表达它的解。分来表达它的解。例例 例例 的解的解的值仍需插值方法来计算的值仍需插值方法来计

5、算.78事实上，从事实上，从实际问题实际问题当中抽象出来的微分方程，通当中抽象出来的微分方程，通常主要依靠常主要依靠数值解法数值解法来解决。来解决。可以证明可以证明:如果函数在带形区域如果函数在带形区域 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上的数值解法上的数值解法。主要讨论一阶

6、常微分方程初值问题主要讨论一阶常微分方程初值问题 9微分方程数值方法的基本思想微分方程数值方法的基本思想 对常微分方程初值问题的数值解法，就是要算对常微分方程初值问题的数值解法，就是要算出精确解出精确解y(x)y(x)在区间在区间 a,ba,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值相邻两个节点的间距相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。假定长可以相等，也可以不等。假定h h为定数，称为定步长，为定数，称为定步长，这时节点可表示为这时节点可表示为10数值解法需要把连续性的问题加以离散化，从数值解法需要把连续性的问题加以离散化，从而求出离散节点的

7、数值解。而求出离散节点的数值解。离散化离散化 对常微分方程数值解法的基本出发点就是对常微分方程数值解法的基本出发点就是离散化离散化。其数值解法的基本特点：采用其数值解法的基本特点：采用“步进式步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，即求解过程顺着节点排列的次序一步一步地向前推进，11描述这类算法，要求给出用已知信息描述这类算法，要求给出用已知信息 计算计算 的递推公式。的递推公式。建立这类递推公式的基本方法是在这些节点上用建立这类递推公式的基本方法是在这些节点上用数数值积分、数值微分、泰勒展开值积分、数值微分、泰勒展开等离散化方法，对初等离散化方法，对初值问题值问题中的导数中

8、的导数 进行不同的进行不同的离散化处理离散化处理。12数值解和精确解数值解和精确解用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。用y(x)表示问题的准确解y(x0),y(x1),y(xN)表示解y(x)在节点x0,x1,xN处的准确值y0,y1,y N表示数值解，即问题的解y(x)在相应节点处的近似值。13单步法和多步法单步法和多步法单步法：在计算yi+1 时只利用y i多步法：在计算yi+1 时不仅利用y i,还要利用 yi1,yi2,k步法：在计算yi+1 时要用到yi,yi1,yik+1显式格式可写成:yk+1=yk+hf(xk，yk;

9、h）隐式格式:yk+1=yk+hf（xk，yk,yk+1;h）它每步求解yk+1需要解一个隐式方程。14欧拉（欧拉（EulerEuler）方法）方法在x=x0 处，用差商代替导数：由由得得15同理，在同理，在x=xx=xn n 处，用差商代替导数：处，用差商代替导数：由由得得若记若记则上式可记为则上式可记为此此即即为为求求解解初初值值问问题题的的EulerEuler方方法法，又又称称显显式式EulerEuler方法。方法。16例例:用用EulerEuler方法求解常微分方程初值问题方法求解常微分方程初值问题并将数值解和该问题的解析解比较。并将数值解和该问题的解析解比较。解：解：EulerEul

10、er方法的具体格式：方法的具体格式：17h=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)程序实现程序实现18xn y(xn)y

11、n 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.计算结果如下：计算结果如下：19xn y(xn)

12、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方法的精度是比较差的。方法的精

13、度是比较差的。20二阶二阶Runge-KuttaRunge-Kutta方法方法其中其中其中其中 c c1 1,c c2 2,2 2,21 21 待定。待定。待定。待定。上式的上式的局部截断误差局部截断误差为：为：21即即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例例 蛇形曲线蛇形曲线的初值问题的初值问题令令f(x,y)=y/x 2yf(x,y)=y/x 2y2 2,取取 f(0,0)=1,f(0,0)=

14、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 23x0=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);f

15、or 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常用的一个公式为常用的一个公式为四阶四阶Runge-KuttaRunge-Kutta方法方法25 function ydot=harmonic(t,y)ydot=y(2);-y(1)y=inline(0 1;-1 0*y,t,y);System of Equations26 function ydot=twobody(t,y)r=sqrt(

16、y(1)2+y(2)2);ydot=y(3);y(4);-y(1)/r3;-y(2)/r3;Two Body Problem27Linearized Differential Equations28J的特征值是解增长解衰减解振荡29 基于龙格库塔法，基于龙格库塔法，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

17、是定义是定义f(t,y)f(t,y)的函数文件名的函数文件名，该函数文件，该函数文件必须返回一个列向量。必须返回一个列向量。tspantspan形式为形式为t0,tf,t0,tf,表示表示求解区间。求解区间。y0y0是初始状态列向量。是初始状态列向量。t t和和y y分别给出分别给出时间向量和相应的状态向量。时间向量和相应的状态向量。MATLABMATLAB求常微分方程数值解的函数求常微分方程数值解的函数30 ode23:Bogacki,Shampine(1989)和Shampine(1994),”23”表示用两算法：一个2阶，一个3阶Bogacki,P.and LF Shampine,A 3

18、(2)pair of Runge-Kutta formulas,Appl.Math.Letters,Vol.2,1989,pp 1-9.BS23 algorithm31 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 ode23(twobody,0 2*pi,1;0;0;1);01234567-1.5-1-0.500.511.5Examples33 y0=1;0;0;3;ode23(twobody,0 2*pi,y0);01234567-2

19、02468101214161834 y0=1;0;0;3;t,y=ode23(twobody,0 2*pi,y0);plot(y(:,1),y(:,2);axis equal35 y0=1;0;0;3;t,y=ode23(twobody,0 2*pi,y0);plot(y(:,1)plot(y(:,2)36 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 t

20、ake small steps to obtain satisfactory results.A model of flame propagation（火焰燃烧）:y是球的半径，y2和y3与球的表面积和体积有关想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小（由于进入球内氧气和消耗的氧气平衡）Stiff Problem(刚性问题)37010203040506070809010000.20.40.60.811.21.4eta=0.02;sym y;F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);38eta=0.00002;sym y;F=

21、inline(y2-y3,t,y);ode23(F,0 2/eta,eta);39012345678910 x 10400.20.40.60.811.21.4ode23seta=0.00002;ode23s(inline(y2-y3,t,y),0 2/eta,eta);40例例 蛇形曲线的常微分方程初值问题蛇形曲线的常微分方程初值问题 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)输出结果为图形输出结果为图形 41T,y=ode23

22、(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

23、 0.5000 0.4868 0.4607 0.4242 0.3793 0.3270 0.2783 0.2411 0.2122 0.1891 0.1705 0.162042例例 洛伦兹模型洛伦兹模型由如下常微分方程组描述由如下常微分方程组描述 取取 =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平

24、面的投影曲线平面的投影曲线 43气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应蝴蝶效应。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前

25、，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。44天气预报的准确性：Lorenz现象的数学：http:/=225分形艺术电子版：混沌理论：http:/www.lamost.org/yzhao/doc/chaos.html45记向量记向量 y y1 1，y y2 2，y y3 3=x x，y y，z z，创建创建MATLABMATLAB函函数文件如下数文件

26、如下 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,

27、y=ode23(flo,0,80,y0)plot(y(:,2),y(:,1)plot(y(:,2),y(:,1)46plot(y(:,3),y(:,1)plot(y(:,3),y(:,1)47plot3(y(:,1),y(:,2),y(:,3)48 y0=30;0;-40;plot(y(:,i)y0=30;0;-40;plot(y(:,i)495051非刚性系统：ode45(Runge-Kutta45)ode23(Runge-Kuatta23)ode113(Adams-Bashforth-Moulton PECE)多步方法刚性系统：ode15s(Gear方法),多步方法 ode23s(二阶mo

28、dified Rosenbroack formula)，单步 ode23t(trapezoidal rule),solve DAEs ode23tb(TR-BDF2)low order methodMatlabs ODE Solvers52Laplacian 算子：Poisson方程（elliptic）：Laplacian 算子的特征值问题:Heat equation(parabolic):Wave equation(hyperbolic):PDE Model53五点离散 Poisson方程离散：特征值问题:Finete Difference Methods54热方程：波动方程：55-4 1

29、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

30、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椭圆方程：特征值方程：热方程：波动方程：离散方程57 例 用古典显式格式求解抛物型方程初始条件为边界条件为取步长x=h=0.2,t=k=0.02。58 解 r=k/h2

31、=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.0259function 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%赋初始条件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%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后再画图626364PDETOOLPDETOOL有限元方法有限元方法65