第14章计算机数学基础下常微分方程的数值解法.doc

第14章 常微分方程的数值解法常微分方程的数值解法的基本思想：求解如下形式的一阶常微分方程： 只有少数具有较简单形式的微分方程才能求出其精确的解析解（解析表达式表示的解）。在实际问题中往往不能求出常微分方程精确的解析解，因此只能利用数值方法来求出微分方程的解在某些点上的近似值，通常称为数值解法。在上述初值问题解的存在区间a,b内，求它在一系列节点：a=x0x1x2x3xn=b上的近似值yk，即yky(xk)(k=1,2,3,n),h=xk+1-xk(k=0,1,2,n-1)称为步长，一般取h为常数(取等距节点)。一、 欧拉法和改进欧拉法1欧拉公式设一阶常微分方程的初值问题：求微分方程根的几何意义：寻找一条经过点（x0,y0）的平面曲线y=y(x)，这条曲线在点(x,y)处的切线斜率是f(x,y)。欧拉折线法的几何意义：用一系列折线来逼近所求的曲线。xy作法：过点(x0,y0)以f(x0,y0)为斜率作切线，其方程为： 当x=x1时，切线的值记为y1,以y1作为y(x1)的近似值,即y(x1)y1,再过点(x1,y1)以f(x1,y1)为斜率作切线，其方程为：重复以上做法，可得到切线的公式为：当x=xn+1时，得到y(xn+1)的一个近似值当所取的节点xk为等距时，有:步长为h=xk+1-xk欧拉公式： 使用欧拉公式求解一阶常微分方程初值问题近似解的方法称为欧拉法（或欧拉折线法）。解题步骤：先明确所给的条件中的已知点(x0,y0)，步长h，以及所给出的f(x,y)，再套用欧拉公式中的递推公式，依次求出y1,y2,y3。例1：取步长h=0.1, 用欧拉法求解初值问题的计算公式解：已知x0=1,y1=1;h=0.1;xk=1+0.1k欧拉公式 此处，迭代公式为：例2：用欧拉法求初值问题在x=0.1,0.2,0.3,1.0处解的近似值。解：已知f(x,y)= , x0=0 , y0=1 , 取等步长h=0.1,n=10由欧拉公式得=yk+0.1（）xk=x0+kh=0.1k（k=0,1,2,9）依次求出解在各点的近似值，把计算过程列表（数值计算表）。欧拉法的局部截断误差：对应于节点xk+1处的精确解为y(xk+1),yk+1是用欧拉公式得到的近似解，则称y(xk+1)yk+1为欧拉公式的局部截断误差。欧拉公式的局部截断误差是O(h2)，表示与h2是同阶无穷小量，即称欧拉法具有一阶精度。例3：用欧拉法解初值问题，取步长h=0.2.解：h=0.2, =yxy2，n=3,xk=0.2k;首先建立欧拉迭代格式： 当k=0，x1=0.2时，已知x0=0,y0=1，有y(0.2)y1=0.21(401)0.8当k1，x2=0.4时，已知x1=0.2, y1=0.8，有y(0.4)y2=0.20.8(40.20.8)0.614 4当k=2，x3=0.6时，已知x2=0.4,y2=0.6144，有y(0.6)y3=0.20.6144(40.40.4613)0.82改进欧拉公式在初值问题中，对方程y=f(x,y)在区间xk,xk+1上直接积分，得到：对上式可用积分梯形公式，即得：用yk+1,yk替代y(xk+1),y(xk)得到计算公式：上式称为求y(xk+1)的近似值的梯形公式。梯形公式并不能直接求得yk+1,称为隐式形式，需要解方程才能得到。因此一般先用欧拉公式求出初始近似值，然后再按迭代方式求解。取等距节点，即h为常数，有也可表示成下式： 或表示成平均形式：上式称为改进的欧拉公式。用改进欧拉公式去求解一价常微分方程初值问题的近似解的方法称为改进欧拉法。改进欧拉公式的局部截断误差：即具有二阶精度，是h3的同阶无穷小量。例4：用欧拉预报校正公式求解初值问题，取步长h=0.2,计算y(0.2),y(0.4)的近似值解： 步长h=0.2, 此时f(x,y)=yy2sinx 欧拉预报校正公式为：有迭代格式：当k=0，x0=1, y0=1时，x1=1.2，有 当k=1，x1=1.2, y1=0.71549时，x2=1.4，有 例5：用改进的欧拉法平均公式，取步长h=0.1，求解初值问题 计算过程保留4位小数。解 首先建立迭代格式:x0=0,y0=1,x1=0.1当k=0时，x0=0,y0=1,x1=0.1,有 当k=1时，x1=0.1, y1=1.11, x2=0.2,有所求y(0.1)1.11; y(0.2)1.242 1 例6：对初值问题，证明：用梯形公式求得的近似解为证明 解初值问题的梯形公式为 整理成显式: 反复迭代，得到 三、龙格库塔法，由微分中值定理可知，存在有0它是一步法；2.精确度高。缺点：计算量大。例7：写出用四阶龙格库塔法求解初值问题的计算公式，取步长h=0.2，计算y(0.4)的近似值，计算过程保留4位小数。解：此处f(x,y)=83y, 四阶龙格库塔法公式为 其中 k1=f(xk,yk)；k2=f(xk+h，yk+hk1)；k3=f(