实验名称微分方程数值解

探索实验8常微分方程初值问题数值解一、实验目的了解常微分方程初值问题的数值解概念，掌握解常微分方程初值问题的Euler方法及改进的Euler方法和Runge-Kutta方法解常微分方程初值问题的算法构造和计算。能用程序实现解常微分方程初值问题的Euler方法、改进的Euler方法和经典的Runge-Kutta方法，学习用计算机求常微分方程初值问题数值解的一些科学计算方法和简单的编程技术。二、概念与结论1.常微分方程初值问题：常微分方程特解问题称为初值问题，通常其形式为2．常微分方程初值问题数值解：常微分方程初值问题的解y(x)在[a,b]上的有限个值y(xk)的近似值yk称为常微分方程初值问题数值解，其中xk=xk-1+hk-1，k=1,2,3,…,N。xk称为节点，hk称为步长。通常，步长h不变，取为等距步长hk=h=（b-a）/N，N为等分区间[a,b]分割数。3.常微分方程初值问题数值解法：求常微分方程初值问题数值解yk的方法称为微分方程初值问题数值解法。4.单步法：在计算yk+1之时只用到yk的方法，其计算公式有：显式单步法计算公式:yk+1=yk+hΦ(xk,yk;h）隐式单步法计算公式:yk+1=yk+hΦ(xk,yk,yk+1,h）式中函数Φ是连续函数，称为增量函数。5.多步法：在计算yk+1之时不仅用到yk，还要用yk-1，yk-2，…,一般m步法要用到yk,yk-1,yk-2，…yk-m+1，多步法也有显式方法和隐式方法之分。6.数值解法的局部截断误差：假设某常微分方程初值问题数值解法在x=xk没有误差，即y(xk)=yk，称Tk+1=y(xk+1)-yk+1为该数值方法的局部截断误差。显式单步法有其局部截断误差为：Tk+1=y(xk+1)-y(xk)-hΦ(xk,y(xk),h）7.数值解法的阶：如果某常微分方程初值问题数值解法的局部截断误差为：Tk+1=(hp+1)则称该数值方法的阶为P,或称该数值方法为P阶方法。数值方法的阶越高，方法越好。三、程序中Mathematica语句解释：1.k++,++k,k--,--k)1()(),(0yaybxayxfyk++表示赋值关系k=k+1,如:k=1;Table[++k,{5}]获得表{2,3,4,5,6}++k表示先处理k的值,再做赋值k=k+1,如:k=1;Table[k++,{5}]获得表{1,2,3,4,5}k--表示赋值关系k=k-1,如:k=1;Table[k--,{5}]获得表{1,0,-1,-2,-3}--k表示先处理k的值,再做赋值k=k-1,如:k=1;Table[--k,{4}]获得表{0,-1,-2,-3}2.x+=k,x*=kx+=k表示x=x+kx*=k表示x=x*k3.For[stat,test,incr,body]以stat为初值,重复计算incr和body直到test为False终止。这里start为初始值,test为条件,incr为循环变量修正式,body为循环体,通常由incr项控制test的变化。注意:上述命令形式中的start可以是由复合表达式提供的多个初值,如果循环体生成Break[]语句,则退出For循环;如果循环体生成Continue[]语句,则由incr的增量进入For语句的下一次循环。四、方法与程序在实际问题中，常常会遇到一些常微分方程初值问题。对这些问题如果只求助于高等数学中介绍的用求通解加确定边界条件的方法去求解往往是无能为力的。因为一方面通解可能求不出来，另一方面所求的解可能是不能用初等函数表达的形式。人们处理这类问题主要采用常微分方程初值问题数值解的方法来求解。这类方法有很多，这里主要介绍解常微分方程初值问题的Euler方法、改进的Euler方法和Runge-Kutta方法的构造过程和算法程序。1.Euler方法Euler方法是最简单的求微分方程数值解的方法。这个方法由于精度不高，实用中较少使用。人们常用Euler方法来说明求微分方程数值解涉及到的一些问题。1.1Euler方法的构造过程：Euler方法涉及的计算公式有很多方法，这里介绍在求微分方程数值解用的较多的Taylor展开法。设f(x,y)充分光滑,xn+1=xn+h,将y(xn+1)在xn点作Taylor展开，得:y(xn+1)=y(xn)+hy(xn)+(h2/2!)y(xn)+(h3)取其关于h的线性部分，有：y(xn+1)y(xn)+hy(xn)注意到y(xn)=f(xn,y(xn)),用yk代替y(xk)并将“”换为等号“=”,得到Euler公式：yn+1=yn+hf(xn,yn)于是，由初始条件y0=y(a),借助Euler公式就可以依次计算出微分方程初值问题的数值解：y1,y2,…,yk称由Euler公式求数值解的方法为Euler方法。显然Euler方法是单步显式方法。因为对Euler公式有局部截断误差为Tn+1=(h2)因此Euler方法是一阶方法。1.2Euler方法算法：1．输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2．计算节点数n和节点xk3．用Euler公式yn+1=yn+hf(xn,yn)求数值解1.3Euler方法程序：Clear[x,y,f]f[x_,y_]=Input[函数f(x,y)=]y0=Input[初值y0=]a=Input[左端点a=]b=Input[右端点b=]h=Input[步长h=]n=(b-a)/h;For[i=0,in,i++,xk=a+i*h;y1=y0+h*f[xk,y0];Print[y(,xk+h//N,)=,y1//N];y0=y1]说明：本程序用Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。程序中变量说明：f[x,y]:存放函数f(x,y)y0:存放初值y0及数值解a：存放自变量区间左端点b:存放自变量区间右端点n:存放节点个数h:存放节点步长xk：存放节点xiy1:存放数值解注：1）语句Print[y(,xk+h//N,)=,y1//N]是将数值解用6位有效数字显示出来，如果要显示n位数有效数字，可以将语句改为Print[y(,xk+h//N,)=,N[y1,n]]2）Mathematica中有求微分方程初值问题数值解的命令，形式为:NDSolve[{y’[x]==f[x,y[x]],y[a]==y0},y,{x,a,b}]由NDSolve命令得到的解是以{{未知函数名-InterpolatingFunction[range,]}}的形式给出的,其中的InterpolatingFunction[range,]是所求的插值函数表示的数值解,range就是所求数值解的自变量范围。1.4例题与实验例1.用Euler方法求初值问题的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。1)0(102yxyxyy解：执行Euler方法程序后，在输入的窗口中按提示分别输入x-2x/y、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=1.1y(0.2)=1.19182y(0.3)=1.27744y(0.4)=1.35821y(0.5)=1.43513y(0.6)=1.50897y(0.7)=1.58034y(0.8)=1.64978y(0.9)=1.71778y(1.)=1.78477本题的准确解为y(x)=(1+2x)1/2，在一个坐标系中画出数值解与准确解的图形如下：图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近，但还是有误差的。2.改进的Euler方法改进的Euler方法比Euler方法精度高，其中把微分方程初值问题数值隐式解法计算公式变为预测、校正公式的做法很有代表性。2.1改进的Euler方法构造过程：将微分方程y=f(x,y)在[xk,xk+1]积分，得对右端的定积分用数值积分的梯形公式，有用yk代替y(xk)并将“”换为等号“=”可得求微分方程初值问题的数值解的梯形公式：11kkxxkkdxxyxfxyxy))(,()()()))(,())(,((2))(,()()(1111kkkkxxkkxyxfxyxfhdxxyxfxyxykk这个公式是单步隐式公式。可以证明它是二阶方法，比Euler方法阶高。不过，在给定初始条件y0=y(a)后,要求出数值解，每一次计算yk的值都要进行非线性方程求根的迭代解法来完成，因此计算量很大。为减少计算量，通常采用迭代一两次就转入下一步计算的方法，特别，如果先用Euler公式计算一次以对要计算的下一步值解进行预测，然后再用梯形公式对其进行校正的方法得到下一步的值，它的计算格式为：这个计算格式称为改进的Euler公式，而称由如上计算格式求数值解的方法为改进的Euler方法。2.2改进的Euler方法算法：1.输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2.计算节点数n和节点xk3.用改进的Euler公式求数值解2.3改进的Euler方法程序：Clear[x,y,f]f[x_,y_]=Input[函数f(x,y)=]y0=Input[初值y0=]a=Input[左端点a=]b=Input[右端点b=]h=Input[步长h=]n=(b-a)/h;For[i=0,in,i++,xk=a+i*h;y1=y0+h*f[xk,y0];y1=y0+h/2*(f[xk,y0]+f[xk+h,y1]);Print[y(,xk+h//N,)=,y1//N];y0=y1]注：语句h=Input[步长h=]的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。说明：本程序用改进的Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。程序中变量说明：f[x,y]:存放函数f(x,y)y0:存放初值y0及数值解a：存放自变量区间左端点b:存放自变量区间右端点n:存放节点个数h:存放节点步长xk：存放节点xi)),(),((1112kkkkkkyxfyxfhyy))ˆ,(),((2),(ˆ1111kkkkkkkkkkyxfyxfhyyyxhfyyy1:存放数值解2.5例题与实验例2.用改进的Euler方法求初值问题的数值解。取步长h=0.1和0.01计算，并与准确解比较。解：执行改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*y^2、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=0.99y(0.2)=0.961366y(0.3)=0.917246y(0.4)=0.861954y(0.5)=0.800034y(0.6)=0.735527y(0.7)=0.671587y(0.8)=0.610399y(0.9)=0.553289y(1.)=0.500919再执行一次改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*y^2、1、0、1、0.05，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.05)=0.9975y(0.1)=0.990087y(0.15)=0.977978y(0.2)=0.961519y(0.25)=0.941158y(0.3)=0.917417y(0.35)=0.890863y(0.4)=0.862076y(0.45)=0.831624y(0.5)=0.800043y(0.55)=0.767819y(0.6)=0.735383y(0.65)=0