第7章-常微分方程的数值解法

第7章常微分方程的数值解法,1引言2Euler法和改进的Euler法3Runge-Kutta法4程序示例,内容提纲（Outline),本章目的要求:,熟练掌握Euler公式，Runge-Kutta公式以及会用以上公式求常微分方程初值问题的数值解与程序设计.,1引言,在工程和科学计算中，所建立的各种常微分方程的初值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。,只含有一个未知数的微分或导数的方程称为常微分方程（OrdinaryDifferentialEquation），形如,考虑初值问题,在高等数学中用解析法求解常微分方程的初值问题,如,又由得,初值问题解为,很多时候解析解求不出来,如,由得,常微分方程数值解的存在唯一性,考虑一阶常微分方程的初值问题,只要f(x,y)在a,bR1上连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L使对任意定义在a,b上的y1(x)和y2(x)都成立，则上述问题存在唯一解。,绝大多数的方程是很难甚至不可能给出解析解的，对于问题(1),要求它的数值解.,2Euler法和改进的Euler法,2.1Euler法(数值微分(折线)法)若将函数y(x)在点xi处的导数y(xi)用差商来表示，即,再用yi近似地代替y(xi)，则初值问题(1)就化为,上式就是所求的Euler公式。,h为步长,一般取等距步长.,2.1Euler法(数值积分法),设将方程的两端从到求积分，即得显然，只要能近似的算出其中的积分项，我们就可以得到计算的公式。若我们用梯形法计算积分项：即可得如下计算公式,的Euler公式为:,h为步长,一般取等距步长.,例1用Euler法求初值问题,的数值解(取h=0.1)。解因为,故由Euler公式得,计算结果:,Euler算法:,h,Euler法求解例1的Matlab程序:,functionEuler(a,b,y0,n)%y0为y的初值n为区间等分数h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);forj=2:n+1y(j)=y(j-1)+h*f(x(j-1),y(j-1);endfork=1:n+1fprintf(x%d=%fty%d=%fn,k-1,x(k),k-1,y(k);endfunctionz=f(x,y)%f为y的一阶导数的函数表达式z=x-y*y;,在编辑窗口输入:,在命令窗口执行:,a=0;b=1;y0=0;n=10;Euler(a,b,y0,n),Euler法的运行结果:,x0=0.000000y0=0.000000 x1=0.100000y1=0.000000 x2=0.200000y2=0.010000 x3=0.300000y3=0.029990 x4=0.400000y4=0.059900 x5=0.500000y5=0.099541x6=0.600000y6=0.148550 x7=0.700000y7=0.206344x8=0.800000y8=0.272086x9=0.900000y9=0.344683x10=1.000000y10=0.422802,2.2改进的Euler法Euler法虽然形式简单，计算方便，但比较粗糙，精度也低。为了达到较高精度的计算公式，对Euler法进行改进.,Euler法:,计算:,用上式替换Euler法,即得改进的Euler法.,改进的Euler公式,这样建立的预报校正系统称为改进的Euler公式：预报校正为方便编程,写成形式：实践表明，改进的Euler公式明显改善了精度.,改进的Euler公式为:,h为步长,一般取等距步长.,改进的Euler算法:,a,用Euler法和改进Euler法解初值问题,解:Euler公式的具体形式为,其中xn=nh=0.1n(n=0,1,10),已知y0=1，由此式可得,例1,取h=0.1,改进的Euler公式进行计算:,x0=0,y0=1时:yp=1+0.1(1-0)=1.1yc=1+0.1(1.1-(2*0.1)/1.1)=1.091818y1=(1.1+1.091818)/2=1.095909,表,改进的Euler算法:,a,用改进Euler法解初值问题,例1,取h=0.1,改进的Euler法求解演示,改进的Euler法求解的Matlab程序:,functionGJeuler(a,b,y0,n)%y0为y的初值n为区间等分数h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);forj=2:n+1yp=y(j-1)+h*f(x(j-1),y(j-1);yc=y(j-1)+h*f(x(j),yp);y(j)=1/2*(yp+yc);endfork=1:n+1fprintf(x%d=%fty%d=%fn,k-1,x(k),k-1,y(k);endfunctionz=f(x,y)%f为y的一阶导数的函数表达式z=y-2*x/y;,在编辑窗口输入:,在命令窗口执行:,a=0;b=1;y0=1;n=10;GJeuler(a,b,y0,n),改进的Euler法运行结果:,x0=0.000000y0=1.000000 x1=0.100000y1=1.095909x2=0.200000y2=1.184097x3=0.300000y3=1.266201x4=0.400000y4=1.343360 x5=0.500000y5=1.416402x6=0.600000y6=1.485956x7=0.700000y7=1.552514x8=0.800000y8=1.616475x9=0.900000y9=1.678166x10=1.000000y10=1.737867,小结:,1.Euler公式,2.改进的Euler公式,3.应用,Runge-Kutta法的设计思想,考察差商，根据微分中值定理，存在点，利用所给方程得我们称为区间上的平均斜率，这样只要对平均斜率提供一种算法，相应地我们便导出一种计算格式。Runge-Kutta方法设计思想就是设法在内多预报几个点的斜率值，然后把它们加权平均作为平均斜率，以期望构造出更高精度的计算格式。,3Runge-Kutta法,二阶Runge-Kutta方法,随意考察区间内一点,用两个点的斜率的加权平均代替平均斜率，于是我们就得到如下计算格式：其中有两个待定参数，适当选取它们的值，就可使上述格式有较高的精度。若，该格式是二阶的，故统称满足这一条件的一族格式为二阶Runge-Kutta格式。特别地，当时，上述格式即为改进的Euler格式，如果取，则上述格式称为变形的Euler格式，亦称为中点格式。,三阶Runge-Kutta方法,为了进一步提高精度，我们可以考虑用三个点的斜率值加权平均得出平均斜率的近似值，其中，于是就可以构造所谓的三阶Runge-Kutta格式，下列Kutta格式是其中的一种：,四阶Runge-Kutta方法,继续上述过程，我们可以导出四阶Runge-Kutta格式，下列经典格式是其中的一种：值得注意的是，Runge-Kutta法的推导基于泰勒展开法，因而它要求解具有较好的光滑性。如果解的光滑性差，则该方法得到的解反而不好。,四、四阶Runge-Kutta方法,例.使用高阶R-K方法计算初值问题,解:,(1)使用三阶R-K方法,其余结果如下:,(2)如果使用四阶R-K方法,nxnk1k2k3yn1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.30001.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993,其余结果如下:,nxnk1k2k3k4yn1.00000.10001.00001.10251.11331.23511.11112.00000.20001.23461.37561.39211.56331.25003.00000.30001.56251.76391.79082.04231.42864.00000.40002.04082.34282.38922.78051.66675.00000.50002.77773.26003.34764.00572.0000,例用四阶RungeKutta法解初值问题,y=x2y(0x1)y(0)=1,解：取h=0.1，,代入计算出相应的数值解即可.,四阶Runge-Kutta方法的Matlab程序:,图6.5,例用四阶RungeKutta法解初值问题的Matlab程序:,y=x2y(0x1)y(0)=1,取h=0.1，,四阶Runge-Kutta方法的Matlab程序:,functionodeRK4(a,b,y0,n)%y0为y的初值n为区间等分数h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);forj=2:n+1k1=f(x(j-1),y(j-1);k2=f(x(j-1)+h/2,y(j-1)+k1*h/2);k3=f(x(j-1)+h/2,y(j-1)+k2*h/2);k4=f(x(j-1)+h,y(j-1)+k3*h);y(j)=y(j-1)+h/6*(k1+k4)+h/3*(k2+k3);endfork=1:n+1fprintf(x%d=%fty%d=%fn,k-1,x(k),k-1,y(k);endfunctionz=f(x,y)%f为y的一阶导数的函数表达式z=x*x-y;,在命令窗口执行:,a=0;b=1;y0=1;n=10;odeRK4(a,b,y0,n),四阶Runge-Kutta方法运行结果:,x0=0.000000y0=1.000000 x1=0.100000y1=0.905163x2=0.200000y2=0.821269x3=0.300000y3=0.749182x4=0.400000y4=0.689680 x5=0.500000y5=0.643470 x6=0.600000y6=0.611189x7=0.700000y7=0.593415x8=0.800000y8=0.590672x9=0.900000y9=0.603431x10=1.000000y10=0.632122,内容大纲,1.常微分方程的数值解(调用内置函数);2.常微分方程的符号解(调用内置函数).,计算实例及其程序演示(调用内置函数),MATLAB软件实现,解析解,dsolve(eqn1,eqn2,c1,var1,),例,输入：y=dsolve(Dy=1+y2)y1=dsolve(Dy=1+y2,y(0)=1,x),输出：y=tan(t-C1)（通解，一簇曲线）y1=tan(x+1/4*pi)（特解，一条曲线）,MATLAB软件实现,例常系数的二阶微分方程,y=dsolve(D2y-2*Dy-3*y=0,x)y=dsolve(D2y-2*Dy-3*y=0,y(0)=1,Dy(0)=0,x),输入：,x=dsolve(D2x-(1-x2)*Dx+x=0,x(0)=3,Dx(0)=0),上述两例的计算结果怎样？由此得出什么结论？,例无解析表达式！,x=dsolve(Dx)2+x2=1,x(0)=0),例非线性微分方程,x=sin(t)-sin(t)若欲求解的某个数值解，如何求解？,t=pi/2;eval(x),MATLAB软件实现,输入：x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y)x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y,x(0)=0,y(0)=1),例,MATLAB软件实现,返回,t，x=solver（f,ts,x0,options）,Matlab软件计算数值解,1）首先建立M-文件（weif.m）functionf=weif(x,y)f=-y+x+1;2）求解：x,y=ode23(weif,0,1,1)3)作图形:plot(x,y,r);4)与精确解进行比较holdonezplot(x+exp(-x),0,1),例1y=-y+x+1，y(0)=1,标准形式：y=f(x,y),范例,1、在解n个未知函数的方程组时，x0和x均为n维向量，m-函数文件中的待解方程组应以x的分量形式写成.,2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.,注意:,注意1：,1、建立M文件函数functionxdot=fun(t,x,y)xdot=f1(t,x(t),y(t)；f2(t,x(t),y(t)；2、数值计算（执行以下命令）t,x,y=ode23(fun,t0,tf,x0,y0),注意：执行命令不能写在M函数文件中。,例如：,令,注意2：,functionxdot=fun1(t,x,y)（fun1.m)xdot=f(t,x(t),y(t)；x(t)；t,x,y=ode23(fun1,t0,tf,x0,y0),M-文件函数如何写呢？,注意：y(t)是原方程的解。x(t)只是中间变量。,如果方程形式是：z=f(t,z,z)?,返回,该方程是否有解析解？,范例,（1）编写M文件(文件名为vdpol.m):functionyp=vdpol(t,y);yp=y(2);(1-y(1)2)*y(2)-y(1);,（2）编写程序如下:（vdj.m）t,y=ode23(vdpol,0,20,3,0);y1=y(:,1);%原方程的解y2=y(:,2);plot(t,y1,t,y2,-)%y1(t),y2(t)曲线图pause,plot(y1,y2),grid,%相轨迹图，即y2(y1)曲线,范例,蓝色曲线y（1）；（原方程解）红色曲线y（2）；,计算结果,范例,范例,例3考虑Lorenz模型：,其中参数=8/3，=10，=28,解：1）编写M函数文件(lorenz.m);2)数值求解并画三维空间的相平面轨线；(ltest.m),范例,1、lorenz.mfunctionxdot=lorenz(t,x)xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1*x;,2、ltest.mx0=000.1;t,x=ode45(lorenz,0,10,x0);plot(t,x(:,1),-,t,x(:,2),*,t,x(:,3),+)pauseplot3(x(:,1),x(:,2),x(:,3),gridon计算结果如下图,范例,图中，x1的图形为实线(蓝），x2的图形为“*”线（绿），x3的图形为“+”线（红）.取t0，tf=0，10。,若自变量区间取0，20、0，40，计算结果如下：,范例,曲线呈震荡发散状,三维图形的混沌状,ltest.m,观察结果：,1、该曲线包含两个“圆盘”，每一个都是由螺线形轨道构成。某些轨道