第8章常微分方程數(shù)值解法

第8章常微分方程數(shù)值解法8.1引言

在常微分方程中利用典型方程的解法，對這些方程就可以求得其解的解析表達式。例如，對一階線性微分方程

(8.1)有通解為有了解析解，就可以根據(jù)初始條件或邊界條件把其中的任意常數(shù)具體確定下來。

在許多工程實踐和科學研究中，除了少數(shù)特殊類型的微分方程能用解析法求得其精確解外，大多數(shù)情況要得出解的解析表達式是極其困難的，甚至是不可能的，因此，只能用近似方法求得其近似解。例如我們很容易求出初值問題的解為但要計算它的值，還需要用數(shù)值積分的方法。如果要對許多個值計算解的近似值，那么工作量非常大。況且實際計算不一定要求解析表達式，而是只需求在某些點上滿足精度的解的近似值或解的近似表達式就可以了。由于高階常微分方程可以轉(zhuǎn)化為一階常微分方程組，因此，為了不失一般性，本章主要介紹一類一階常微分方程初值問題

在區(qū)間上的數(shù)值計算方法。在初值問題的數(shù)值解中，常見的有單步法和多步法。所謂單步法，就是在算法中根據(jù)初始條件，由解就能求出；所謂多步法，就是在求時不僅要知道，而且還要知道在前面的值。本章首先介紹單步法的歐拉法及龍格-庫塔法，然后介紹多步法的阿當姆斯法，它們是解初值問題的幾個最基本的常用的有效方法。8.2歐拉（Euler）法

8.2.1歐拉法的基本思想雖然歐拉法在實際應用時較少采用，但由于它的結(jié)構(gòu)簡單，在理論上具有非常重要的意義。歐拉法的基本思想就是用差分方程初值問題(8.3)的解來近似微分方程初值問題(8.2)的解，其中，式(8.3)也稱為歐拉公式。

歐拉法的幾何意義是用一條自點出發(fā)的折線去逼近積分曲線，如圖8.1所示。因此，這種方法又稱為折線法。顯然，歐拉法簡單地取折線的端點作為數(shù)值解，精度非常差。

8.2.2實現(xiàn)歐拉法的基本步驟（1）輸入的區(qū)間，區(qū)間的等分數(shù)以及在處的初始值；（2）對，計算：

（3）輸出；（4）結(jié)束。例8.1用歐拉法求解初值問題的數(shù)值解（保留4位小數(shù)）。#include<stdio.h>floatfunc(floatx,floaty){return(y-x);}floateuler(floatx0,floatxn,floaty0,intN){floatx,y,h;

inti;x=x0;y=y0;h=(xn-x0)/(float)N;/*計算步長*/

for(i=1;i<=N;i++)/*歐拉公式*/{y=y+h*func(x,y);x=x0+i*h;}return(y);}main(){floatx0,xn,y0,e;

intN;

clrscr();

printf("\ninputn:\n");/*輸入?yún)^(qū)間等分數(shù)*/

scanf("%d",&N);

printf("inputx0,xn:\n");/*輸入x的區(qū)間*/

scanf("%f%f",&x0,&xn);

printf("inputy0:\n");/*輸入x0處的y的值*/

scanf("%f",&y0);e=euler(x0,xn,y0,N);

printf("y(%f)=%6.4f",y0,e);}程序運行結(jié)果：inputn:5inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.4883inputn:10inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.5937inputn:20inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.6533此問題的精確解為，相應計算值為4.7183。從計算結(jié)果可以看出，雖然等分數(shù)越大，數(shù)值解越精確，但精度并不理想。8.3改進的歐拉法

8.3.1改進的歐拉法的基本思想改進的歐拉法以隱函數(shù)的方式給出，計算時常使用迭代法，一般可先由歐拉公式(8.3)計算出的初始值，再按下面的格式進行的迭代計算(8.4)改進的歐拉法由于添加了迭代過程，計算量較歐拉法增加了一倍，付出這種代價的目的是為了提高精度。8.3.2實現(xiàn)改進的歐拉法的基本步驟（1）輸入的區(qū)間，區(qū)間的等分數(shù)以及在處的初始值；（2）,，；（3）；（4）計算：（5）若，，轉(zhuǎn)（4）；否則，轉(zhuǎn)（6）；（6）輸出；（7）結(jié)束。例8.2用改進的歐拉法求解初值問題的數(shù)值

解（保留4位小數(shù)）。#include<stdio.h>floatfunc(floatx,floaty){return(y-x);}floateuler(floatx0,floatxn,floaty0,intN){floatx,y,yp,yc,h;

inti;x=x0;y=y0;h=(xn-x0)/(float)N;/*計算步長*/

for(i=1;i<=N;i++)/*改進的歐拉公式*/{yp=y+h*func(x,y);x=x0+i*h;

yc=y+h*func(x,yp);y=(yp+yc)/2.0;}return(y);}main(){floatx0,xn,y0,e;

intN;

clrscr();

printf("\ninputn:\n");/*輸入?yún)^(qū)間等分數(shù)*/

scanf("%d",&N);

printf("inputx0,xn:\n");/*輸入x的區(qū)間*/

scanf("%f%f",&x0,&xn);

printf("inputy0:\n");/*輸入x0處的y的值*/

scanf("%f",&y0);e=euler(x0,xn,y0,N);

printf("y(%f)=%6.4f",y0,e);}程序運行結(jié)果：inputn:5inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.7027inputn:10inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.7141inputn:20inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.71728.4龍格—庫塔（Runge-Kutta）法

8.4.1龍格—庫塔法的基本思想在歐拉法中，用解函數(shù)在點處的斜率計算從到的增量，的表達式與的Taylor展開式的前二項相等，使方法只有一階精度。改進的歐拉法用兩個點，處的斜率、的平均值計算增量，使方法具有二階精度，即的表達式與的Taylor展開式的前三項相等。由此龍格和庫塔提出了一種間接地運用Taylor公式的方法，即利用在若干個待定點上的函數(shù)值和導數(shù)值做出線性組合式，選取適當系數(shù)使這個組合式進行Taylor展開后與的Taylor展開式有較多的項達到一致，從而得出較高階的數(shù)值公式，這就是龍格—庫塔法的基本思想。8.4.2實現(xiàn)標準四階龍格—庫塔法的基本步驟（1）輸入的區(qū)間，區(qū)間的等分數(shù)以及在處的初始值；（2）,，；（3）；（4）計算：（5）若，，轉(zhuǎn)（4）；否則，轉(zhuǎn)（6）；（6）輸出；（7）結(jié)束。例8.3用標準四階龍格—庫塔法求解初值問題的數(shù)值解（保留6位小數(shù)）。#include<stdio.h>floatfunc(floatx,floaty){return(x-y);}floatrunge_kutta(floatx0,floatxn,floaty0,intN){floatx,y,y1,y2,h,xh;floatd1,d2,d3,d4;

inti;x=x0;y=y0;h=(xn-x0)/(float)N;/*計算步長*/

for(i=1;i<=N;i++)/*標準四階龍格-庫塔公式*/{xh=x+h/2;d1=func(x,y);d2=func(xh,y+h*d1/2.0);d3=func(xh,y+h*d2/2.0);d4=func(xh,y+h*d3);y=y+h*(d1+2*d2+2*d3+d4)/6.0;x=x0+i*h;}return(y);}

main(){floatx0,xn,y0,e;

intN;

clrscr();

printf("\ninputn:\n");/*輸入?yún)^(qū)間等分數(shù)*/

scanf("%d",&N);

printf("inputx0,xn:\n");/*輸入x的區(qū)間*/

scanf("%f%f",&x0,&xn);

printf("inputy0:\n");/*輸入x0處的y的值*/

scanf("%f",&y0);e=runge_kutta(x0,xn,y0,N);

printf("y(%f)=%8.6f",y0,e);}程序運行結(jié)果：inputn:5inputx0,xn:0.01.0inputy0:0.0y(0.000000)=0.356261inputn:10inputx0,xn:0.01.0inputy0:0.0y(0.000000)=0.362344inputn:20inputx0,xn:0.01.0inputy0:0.0y(0.000000)=0.365179注意：

此問題的精確解為0.367879。從計算結(jié)果可以看出，隨著等分數(shù)的增加，數(shù)值解越精確。但需要指出的是：龍格—庫塔法的推導是在用Taylor展開法的基礎(chǔ)上進行的，因而它要求所求的解具有較好的光滑性質(zhì)。假若解的光滑性差，那么使用標準四階龍格—庫塔法求得的數(shù)值解，其精度可能反而不如改進的歐拉法。因此，在實際計算時，應當針對問題的具體特點，選擇合適的算法。8.5阿當姆斯（Adams）法

8.5.1阿當姆斯法的基本思想

龍格—庫塔法是一類重要算法，但是，這類算法在每一步都要先預報幾個點上的斜率值，因而計算量比較大。事實上，考慮到在計算之前已得出一系列節(jié)點上的斜率值，能否利用這些已知信息來減少計算量，這正是阿當姆斯法的基本思想。構(gòu)成阿當姆斯法的公式有阿當姆斯外插公式和阿當姆斯預測—校正公式兩種，本節(jié)主要介紹阿當姆斯預測—校正公式。

8.5.2實現(xiàn)阿當姆斯法的基本步驟（1）輸入的區(qū)間，區(qū)間的等分數(shù)以及在處的初始值；（2）；（3）使用標準四階龍格—庫塔法計算,為阿當姆斯預測—校正公式提供初始值；（4）（5）；（6）若，，轉(zhuǎn)（5）；否則，轉(zhuǎn)（7）；（7）輸出；（8）結(jié)束。例8.4用阿當姆斯法求解初值問題的數(shù)值解（保留4位小數(shù)）。#include<stdio.h>floatfunc(floatx,floaty){return(y-2.0*x/y);}voidrunge_kutta(floatx0,floaty0,floath,floatyy[3]){floatx,y,y1,y2,xh;floatd1,d2,d3,d4;

inti;

x=x0;y=y0;for(i=1;i<=3;i++){xh=x+h/2;d1=func(x,y);d2=func(xh,y+h*d1/2);d3=func(xh,y+h*d2/2);d4=func(xh,y+h*d3);y=y+h*(d1+2*d2+2*d3+d4)/6;x=x0+i*h;

yy[i-1]=y;}}floatadams(floatx0,floatxn,floaty0,intN){floath,x,y,yy[3];floaty1,y2,yp,yc;

inti;h=(xn-x0)/(float)N;/*計算步長*/

runge_kutta(x0,y0,h,yy);/*標準四階龍格-庫塔公式*/

y1=yy[0];y2=yy[1];y=yy[2];for(i=3;i<N;i++){x=x0+i*h;

yp=y+h*(55*func(x,y)-59*func(x-h,y2)+37*func(x-2*h,y1)-9*func(x-3*h,y0))/24.0;

yc=y+h*(9*func(x+h,yp)+19*func(x,y)-5*func(x-h,y2)+func(x-2*h,y1))/24.0;y0=y1;y1=y2;y2=y;y=yc;}return(yc);}

main(){floatx0,xn,y0,ad;

intN;

clrscr();

printf("\ninputn:\n");/*輸入?yún)^(qū)間等分數(shù)N*/

scanf("%d",&N);

printf("inputx0,xn:\n");/*輸入x的區(qū)間*/

scanf("%f%f",&x0,&xn);

printf("inputy0:\n");/*輸入x0處的y的值*/

scanf("%f",&y0);ad=adams(x0,xn,y0,N);

printf("y(%f)=%6.4f",y0,ad);}

程序運行結(jié)果：inputn:5inputx0,xn:0.01.0inputy0:1.0y(1.000000)=1.7722inputn:10inputx0,xn:0.01.0inputy0:1.0y(1.000000)=1.7465inputn:20inputx0,xn:0.01.0inputy0:1.0y(1.000000)=1.7364注意：

此問題的解析解為1.7321。從計算結(jié)果可以看出，隨著等分數(shù)的增加，數(shù)值解越精確。阿當姆斯預測—校正公式的計算量比龍格—庫塔法少，但是它的缺點是必須用其它方法（本節(jié)使用標準四階龍格—庫塔法）求出初始值y[0],yy[1],yy[2],因而它不能自啟動。本章小結(jié)

本章主要介紹了四種在計算機上比較常