《解微分方程.ppt》由会员分享，可在线阅读，更多相关《解微分方程.ppt（65页珍藏版）》请在三一办公上搜索。
1、1,MATLABODE初值问题的数值解PDE问题的数值解,2,问题提出 倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式,其中直径D为常数.对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下表是经过测量得到部分容器高度与直径的关系.,3,x,1,o,根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴,问题为如何根据任意高度x标出容器体积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表示高度,直径D是高度x的函数,记为D
2、(x).,4,只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是函数D(x)无解析表达式,无法求出其解析解.,因此,得到如下微分方程初值问题,5,包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。,微分方程分类,6,常微分方程：,(1)，(2)式称为初值问题.,在实际应用中还经常需要求解常微分方程组：,(3)式称为边值问题。,7,但能求解析解的常微分方程是有限的，大多数的常微分
3、方程是给不出解析解的.,这个一阶微分方程就不能用初等函数及其积分来表达它的解。,例,例,的解,的值仍需插值方法来计算.,8,9,事实上，从实际问题当中抽象出来的微分方程，通常主要依靠数值解法来解决。,10,微分方程数值方法的基本思想 对常微分方程初值问题的数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值,相邻两个节点的间距 称为步长，步长可以相等，也可以不等。假定h为定数，称为定步长，这时节点可表示为,11,数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。,离散化,对常微分方程数值解法的基本出发点就是离散化。其数值解法的基本特点：采用“步进式”，即求
4、解过程顺着节点排列的次序一步一步地向前推进，,12,描述这类算法，要求给出用已知信息 计算 的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题,中的导数 进行不同的离散化处理。,13,数值解和精确解,用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。,用y(x)表示问题的准确解,y(x0),y(x1),y(xN)表示解y(x)在节点x0,x1,xN处的准确值y0,y1,y N表示数值解，即问题的解y(x)在相应节点处的近似值。,14,单步法和多步法,单步法：在计算yi+1 时只利用y i多步
5、法：在计算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需要解一个隐式方程。,15,欧拉（Euler）方法,在x=x0 处，用差商代替导数：,由,得,16,同理，在x=xn 处，用差商代替导数：,由,得,若记,则上式可记为,此即为求解初值问题的Euler方法，又称显式Euler方法。,17,例:用Euler方法求解常微分方程初值问题,并将数值解和该问题的解析解比较。,解：Euler方法的具体格式：,1
6、8,h=0.2;y(1)=0.2;x=0.2:h:3;for n=1:14 xn=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),程序实现,19,xn 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.47
7、300.53540.0624,h=0.2,xn=nh,(n=0,1,2,15),f(x,y)=y/x 2y2 计算中取f(0,0)=1.计算结果如下：,20,xn 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,由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小
8、数点后第二位或第三位小数上，说明Euler方法的精度是比较差的。,21,二阶Runge-Kutta方法,其中 c1,c2,2,21 待定。,上式的局部截断误差为：,22,即,c1=1-a，2=21=1/(2a),方程组解不唯一，可令c2=a 0，则,满足上述条件的公式都为2阶R-K公式。,23,例 蛇形曲线的初值问题,令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)2 yn+1=yn+0.5h k1+k2,24,x0=0;y0=0;
9、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);end,y1=x./(1+x.2);plot(x,y,o,x,y1),25,常用的一个公式为,四阶Runge-Kutta方法,26,function ydot=harmonic(t,y)ydot=y(2);-y(1),y=inline(0 1;-1 0*y,t,y)
10、;,System of Equations,27,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 Problem,28,Linearized Differential Equations,29,J的特征值是,解增长,解衰减,解振荡,30,基于龙格库塔法，MATLAB求常微分方程数值解的函数，一般调用格式为：t,y=ode23(fname,tspan,y0)t,y=ode45(fname,tspan,y0)其中fname是定义f(t,y)的函数文件名，该函数文件必须返
11、回一个列向量。tspan形式为t0,tf,表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。,MATLAB求常微分方程数值解的函数,31,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 algorithm,32,F=inline(y(2);-y(1),t,y)ode23(F,0 2*p
12、i,1;0),opts=odeset(reltol,1.e-4,abstol,1.e-6,outputfcn),Examples,33,ode23(twobody,0 2*pi,1;0;0;1);,Examples,34,y0=1;0;0;3;ode23(twobody,0 2*pi,y0);,35,y0=1;0;0;3;t,y=ode23(twobody,0 2*pi,y0);plot(y(:,1),y(:,2);axis equal,36,y0=1;0;0;3;t,y=ode23(twobody,0 2*pi,y0);plot(y(:,1)plot(y(:,2),37,A problem
13、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与球的表面积和体积有关,想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小（由于进入球内氧气和消耗的氧气平衡）,St
14、iff Problem(刚性问题),38,eta=0.02;sym y;F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);,39,eta=0.00002;sym y;F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);,40,eta=0.00002;ode23s(inline(y2-y3,t,y),0 2/eta,eta);,41,例 蛇形曲线的常微分方程初值问题,MATLAB数值求解命令F=inline(1./(1+x.2)-2*y.2)；ode23(F,0,6,0)输出结果为图形,42,T,y=ode23(f,0,6,0)将得到
15、自变量和函数的离散数据,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.0000,y=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
16、0.2783 0.2411 0.2122 0.1891 0.1705 0.1620,43,例 洛伦兹模型由如下常微分方程组描述,取=8/3，=10，=28。初值：x(0)=0，y(0)=0，z(0)=0.01。利用MATLAB求解常微分方程数值解命令计算出t0，80内，三个未知函数的数据值，并绘出相空间在Y-X平面的投影曲线,44,气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平
17、时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。,45,天气预报
18、的准确性：http:/.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLorenz现象的数学：http:/www.sciscape.org/news_detail.php?news_id=225分形艺术电子版：http:/,46,记向量 y1，y2，y3=x，y，z，创建MATLAB函数文件如下,function z=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(
19、3);,用MATLAB命令求解并绘出Y-X平面的投影图,y0=0;0;0.01;x,y=ode23(flo,0,80,y0)plot(y(:,2),y(:,1),47,plot(y(:,3),y(:,1),48,plot3(y(:,1),y(:,2),y(:,3),49,y0=30;0;-40;plot(y(:,i),50,51,52,非刚性系统：,ode45(Runge-Kutta45)ode23(Runge-Kuatta23)ode113(Adams-Bashforth-Moulton PECE)多步方法,刚性系统：,ode15s(Gear方法),多步方法 ode23s(二阶modifie
20、d Rosenbroack formula)，单步 ode23t(trapezoidal rule),solve DAEs ode23tb(TR-BDF2)low order method,Matlabs ODE Solvers,53,Laplacian 算子：,Poisson方程（elliptic）：,Laplacian 算子的特征值问题:,Heat equation(parabolic):,Wave equation(hyperbolic):,PDE Model,54,五点离散,Poisson方程离散：,特征值问题:,Finete Difference Methods,55,热方程：,波动
21、方程：,56,Matrixs,57,椭圆方程：,特征值方程：,热方程：,波动方程：,离散方程,58,例 用古典显式格式求解抛物型方程初始条件为边界条件为取步长x=h=0.2,t=k=0.02。,59,解 r=k/h2=0.02/0.22=0.5,古典显式格式为,0,0 0.2,0 0.4,0 1,0,0,0.20,0.02,1,0.21,0.02,60,function 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
22、y=fg(x)y=4.*(x-x.2);,61,%赋初始条件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);endend,62,%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后再画图,63,64,65,PDETOOL,有限元方法,
