计算方法第8章 常微分方程初值问题的数值解法

4
8.2 欧拉方法
一、显式欧拉法
在一阶常微分方程初值问题的数值解法中，显式欧拉（Euler）法是最简 单的一种。显式欧拉法有明显的几何含义，缺点是精度不高。对于一阶 常微分方程的初值问题：
dy f (x, y) , a x b dx y(x0 ) y0
显式欧拉法的递推公式为：yk=yk-1＋hf(xk-1,yk-1)，k=1,2,3,……。 显式欧拉法每一轮递推只用到前面一轮递推的结果，因此它是单步法。
间[xk-1,xk+1]的中点xk。
12
8.2 欧拉方法
三、改进的欧拉法
改进的欧拉法是一种预测—校正方法，它的每一轮递推包括预测和校正 这2个步骤：先用显式欧拉公式计算出y*k，
y*k=yk-1＋hf(xk-1,yk-1) 这一步称为“预测”；再用梯形公式迭代一次，
yk=yk-1＋(h/2)*(f(xk-1,yk-1)+f(xk,y*k)) 这一步称为“校正”。改进的欧拉法精度比显式欧拉法高，不需要解方 程，是一种更实用的方法。
void main(void)
}
{
double a,b,h,x[MAXSIZE],y[MAXSIZE];
long i,n;
printf("\n请输入求解区间a,b：");
scanf("%lf,%lf",&a,&b);
printf("\n请输入步长h："); scanf("%lf",&h);
n=(long)((b-a)/h);
由基于数值微分的离散化方法、基于数值积分的离散化方法、基于泰 勒展开的离散化方法都可以推导出显式欧拉法的递推公式。
用显式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的(x0,y0) 作为起点，代入显式欧拉法的递推公式的右端，计算出y(x)在x1处的近似 值y1；再以(x1,y1)作为起点，用显式欧拉法的递推公式计算出y(x)在x2处 的近似值y2；……。
x[0]=a;
printf("\n请输入起点x[0]=%lf处的纵坐标y[0]：",x[0]);
scanf("%lf",&y[0]);
for(i=0;i<n;i++)
{
x[i+1]=x[i]+h;
y[i+1]=y[i]+h*f(x[i],y[i]); }
printf("\n计算结果为：");
for(i=0;i<=n;i++)
定义：假设某一步递推的起点是精确的，即yi-1=y(xi-1)，那么这一步递推 的截断误差Ri=y(xi)－yi称为局部截断误差。 定义：若某算法的局部截断误差为O(hp+1)，则称此算法有p阶精度。
显式欧拉法具有一阶精度。
显式欧拉法在步长过大时误差较大；在步长较小时需要多步递推，可能 出现误差积累的现象。由于显式欧拉法的精度不高，因此在实际应用中 用的较少。
dy f (x, y) , a x b dx y(x0 ) y0
目前在常微分方程理论中，只能求出某些特殊类型常微分方程的解析解， 对大部分常微分方程，用解析方法求出常微分方程的精确解非常困难， 甚至不存在解的解析表达式。为满足工程实践的需要，常常用数值解法 求常微分方程的近似解。
在本章中假设讨论的一阶常微分方程的初值问题的解y(x)存在、唯一且 足够光滑，方程本身是稳定的，即精确解y(x)连续且依赖于初始值及右 端函数。
一阶常微分方程初值问题的数值解法有多种分类方法。 一种分类方法为： ⑴ 单步法：每一轮递推只用到前面一轮的递推结果，递推格式为： yk=yk-1＋hT(xk-1,yk-1) ⑵ 多步法：每一轮递推要用到前面多轮递推的结果，递推格式为： yk=yk-1＋hT(xk-r,yk-r,xk-r+1,yk-r+1,……, xk-1,yk-1)，其中r＞1。
13
8.2 欧拉方法
改进的欧拉法的算法
输入求解区间的边界a,b，步长h，起点x[0]=a处的纵坐标y[0]。
递推的次数n=(b-a)/h。
for(i=0;i<n;i++)
//循环1次完成1轮递推
x[i+1]=x[i]+h;
显式欧拉法把ξk近似为区间[xk-1,xk]的起点xk-1； 隐式欧拉法把ξk近似为区间[xk-1,xk]的终点xk； 梯形公式法把处ξk的导数y'(ξk)近似为区间[xk-1,xk]的起点和终点导数
的平均值；
中点欧拉方法考察的区间为[xk-1,xk+1]，ξk∈(xk-1,xk+1)，ξk被近似为区
2
8.1 基础知识
二、数值解法
一阶常微分方程初值问题的数值解法的主要思想，是对区间[a,b]上的节点 a=x0＜x1＜……＜xn＜xn+1＜……≤b
建立y(xn)的近似值yn的某一递推格式，利用初值y0和已计算出的y1,y2,……, yk-1递推出yk，并且用这一方法反复递推，依次得到yk+1,yk+2,……,yn。这一 求解方法称为步进式求解，相邻2个节点的距离称为步长，记为hi=xi+1－xi。 为便于计算，常取成等距节点，称为定步长，这时把步长记为h。
6
8.2 欧拉方法
显式欧拉法的算法
输入求解区间的边界a,b，步长h，起点x[0]=a处的纵坐标y[0]。
递推的次数n=(b-a)/h。
for(i=0;i<n;i++)
//循环1次完成1轮递推
x[i+1]=x[i]+h;
来自百度文库
y[i+1]=y[i]+h*f(x[i],y[i]);
for(i=0;i<=n;i++) 输出x[i],y[i]。
//输出计算结果
7
8.2 欧拉方法 显式欧拉法对应的程序
#include <stdio.h>
double f(double x,double y)
#define MAXSIZE 50 {
double f(double x,double y);
return(……); /*计算并返回函数值f(x,y)*/
10
8.2 欧拉方法
二、欧拉方法的变形
方法三：中点欧拉方法（两步欧拉方法 ） 递推公式为： yk+1=yk-1＋2hf(xk,yk)，k=1,2,3,……。 中点欧拉方法是双步法，需要2个初值y0和y1才能启动递推过程。一般 先用单步法由点(x0,y0)计算出(x1,y1)，再用中点欧拉方法反复地递推。 中点欧拉方法由点(x0,y0)和(x1,y1)递推出点(x2,y2)的几何含义，是经过 点(x0,y0)做一直线段，在区间[x0,x2]内用此直线段代替y(x)的曲线段。 此直线段的斜率等于曲线y(x)在点(x1,y1)处的切线斜率。类似地，由 点(xk-1,yk-1)和点(xk,yk)递推出点(xk+1,yk+1) ，其中k=1,2,3,……。 中点欧拉方法具有2阶精度。
第8章 常微分方程初值问 题的数值解法
❖基础知识 ❖欧拉方法 ❖龙格—库塔方法
1
8.1 基础知识
一、问题的提出
很多实际问题都需要求解常微分方程。例如单摆问题。
常微分方程分为线性常微分方程和非线性常微分方程，又可以分为一阶 常微分方程和高阶常微分方程。通过变量的替换，可以把高阶常微分方 程转化为一阶常微分方程再求解。对于一阶常微分方程组，可以写成向 量形式的单个方程，求解方法与一阶常微分方程相似。因此本章只讨论 一阶常微分方程的初值问题：
⑶ 预测—校正系统：每一轮递推包括预测和校正这2个步骤。先用显式方 法计算出yk，作为迭代的初值，这一过程称为预测；再把隐式方法的递推 公式作为迭代公式，把预测值yk代入迭代公式右端进行迭代，这一过程称 为校正。在校正时往往迭代1次或几次，校正值的精度就会有大幅提高。 一阶常微分方程初值问题的数值解法一般是对连续的初值问题进行离散化 处理，把微分方程转化为代数方程来求解。常用的离散化方法有： ⑴ 基 于数值微分的离散化方法，⑵ 基于数值积分的离散化方法，⑶ 基于泰勒 展开的离散化方法。
printf("\nx[%ld]=%lf,y[%ld]=%lf",i,x[i],i,y[i]);}
8
8.2 欧拉方法
二、欧拉方法的变形
方法一：隐式欧拉法（后退的欧拉法）
递推公式为： yk=yk-1＋hf(xk,yk)，k=1,2,3,……。 用隐式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的 (x0,y0)作为起点，计算出y(x)在x1处的近似值y1；再以(x1,y1)作为起点， 计算出y(x)在x2处的近似值y2；……，象这样反复地递推。 显然，隐式欧拉法是隐式方法，递推公式的右端有未知量yk。隐式欧拉 法需要求解1个方程。为避免解方程，常用显式欧拉法的计算结果作为迭 代的初值yk(0)，把隐式欧拉法的递推公式作为迭代公式反复迭代，得到 迭代序列yk(0),yk(1),yk(2),……。如果步长足够小，那么迭代序列收敛于yk。 隐式欧拉法由点(x0,y0)递推到点(x1,y1)的几何含义，是把曲线y(x)在点 (x1,y1)处的切线平行移动，移动到经过点(x0,y0) ，在区间[x0,x1]内用此 直线段代替y(x)的曲线段。
隐式欧拉法具有一阶精度。
隐式欧拉法精度不高，计算复杂，用的比较少。
9
8.2 欧拉方法
二、欧拉方法的变形
方法二：梯形公式法
递推公式为： yk=yk-1＋(h/2)*(f(xk-1,yk-1)+f(xk,yk))，k=1,2,3,……。 用梯形公式法求解一阶常微分方程初值问题的过程，就是以已知的 (x0,y0)作为起点，计算出y(x)在x1处的近似值y1；再以(x1,y1)作为起点， 计算出y(x)在x2处的近似值y2；……，象这样反复地递推。 显然，梯形公式法是隐式方法，需要求解方程。为避免解方程，常用 显式欧拉法的计算结果作为迭代的初值yk(0)，把梯形公式法的递推公 式作为迭代公式反复迭代，得到迭代序列yk(0),yk(1),yk(2),……。如果步 长h足够小，那么迭代序列收敛于yk。 梯形公式法由点(x0,y0)递推到点(x`1,y1)的几何含义，是经过点(x0,y0)做一 直线段，在区间[x0,x1]内用此直线段代替y(x)的曲线段。此直线段的斜 率等于曲线y(x)在点(x0,y0)处的切线斜率和曲线y(x)在点(x1,y1)处的切线 斜率的平均值。类似地，由点(xk-1,yk-1)递推出点(xk,yk)，k=1,2,3,……。 梯形公式法具有2阶精度。
5
8.2 欧拉方法
一、显式欧拉法（续）
显式欧拉法的递推公式是斜率为f(xk-1,yk-1)，经过点(xk-1,yk-1)的直线方程。 上述显式欧拉法递推过程的几何含义，就是用曲线y(x)在点(x0,y0)处的切 线段代替y(x)在区间[x0,x1]内的曲线段；再把曲线段的终点(x1,y(x1))近似 为切线段的终点(x1,y1)，把曲线y(x)在点(x1,y1)处切线的斜率f(x1,y(x1))近 似为f(x1,y1) ，做曲线y(x)在点(x1,y1)处的切线，并在区间[x1,x2]内用近似 的切线段代替y(x)的曲线段；……。欧拉法用一系列折线近似地代替曲线 y(x)，因此它又称为欧拉折线法。
多步法不能自行启动，必须先用单步法计算出y1,y2,……,yk-1，才能启动一 个r步的多步法。
3
8.1 基础知识
二、数值解法（续）
另一种分类方法为：
⑴ 显式方法：递推公式的右端都是已知量，可以直接计算出递推的结 果，递推格式为：yk=yk-1＋hT(xk-r,yk-r,xk-r+1,yk-r+1,……, xk-1,yk-1) ⑵ 隐式方法：递推公式左端的未知量也出现在公式的右端，递推格式 为： yk=yk-1＋hT(xk-r,yk-r,xk-r+1,yk-r+1,……, xk,yk) 隐式方法的递推公式其实是一个方程。解方程的运算量可能较大，为避免 解方程，常采用预测—校正系统。
11
8.2 欧拉方法
二、欧拉方法的变形
上述各种方法的比较：
由拉格郎日中值定理，在区间(xk-1,xk)内必定存在ξk，使得： y(xk)－ y(xk-1)=y'(ξk)(xk－xk-1)
∴y(xk)=y(xk-1)+hy'(ξk)=y(xk-1)+hf(ξk,y(ξk))，
如果知道ξk，代入这个递推公式，那么递推过程得到的序列y0,y1, y2,……没有误差。求ξk往往很困难，因此常用一个易求的值近似地代替 ξk 。显式欧拉法、隐式欧拉法、梯形公式法、中点欧拉方法的区别是 对ξk 、y'(ξk)的近似方法不同。