工程计算基础第9章-1

1、第10章 常微分方程的数值解法微分方程 微分方程是精确表示自然科学中各种基本定律和各种问题的基本工具之一。现代建立起来的自然科学和社会科学中的数学模型大多都是微分方程。 在许多物理、力学、生物等现象中，不能直接找到联系所研究的那些量的规律，但却容易建立起这些量与它们的导数或微分间的关系。含有未知函数的导数（或微分）的方程，称为微分方程。未知函数可以不出现，但其导数一定要出现。未知函数可以不出现，但其导数一定要出现。车辆纵向动力学模型mFa直流电机模型牛顿定律和基尔霍夫定律电机特性其中： T为电磁转矩； i为电枢电流； Kt为转矩系数； e为反电动势； Ke为反电动势系数。二自由度转向模型（1）

2、忽略转向系统的影响，直接以前轮转角作为输入；（2）忽略悬架的作用；车身只作平行于地面的平面运动，沿z 轴的位移、绕 y轴的俯仰角和绕 x 轴的侧倾角均为零，且lrZZFF（3）汽车前进速度u视为不变；（4）侧向加速度限定在0.4g一下，确保轮胎侧偏特性处于线性范围； （5）驱动力不大，不考虑地面切向力对轮胎侧偏特性的影响，没有空气动力的作用。常微分方程 什么是常微分方程 ODE，Ordinary Differential Equations(1) y= kx, k 为常数；(2) ( y - 2xy) dx + x2 dy = 0；(3) mv(t) = mg - kv(t)；;112yay

3、(4).,(0sind22为常数lglgtd(5)常微分方程、偏微分方程、微分方程的阶 未知函数是一元函数的微分方程称做常微分方程， 未知函数是多元函数的微分方程称做偏微分方程 微分方程中出现的未知函数最高阶导数的阶数，称为微分方程的阶 n 阶微分方程的一般形式为：F(x, y, y, , y(n) = 0，其中 x 是自变量， y 是未知函数，F(x, y, y, , y(n)是已知函数，222222222y=4xy -(y) +12xy=000ydydytydtdtTTTxyz 一阶线性二阶非线性一阶非线性二阶偏微分方程例如例如关于微分方程的解 任何代入微分方程后使其成为恒等式的函数，都叫

4、做该方程的解. 若微分方程的解中含有任意常数的个数与方程的阶数相同，且任意常数之间不能合并，则称此解为该方程的通解(或一般解). 当通解中的各任意常数都取特定值时所得到的解，称为方程的特解. 为了确定特解中的任意常数，对n阶方程，相应必须提出n个条件，这样的条件称为定解条件。微分方程与定解条件联立构成定解问题。微分方程的初值问题和边值问题 给出方程定义区间左端点处的某些值（初始值），这样的条件称为初始条件。方程和初始条件构成的问题称为初值问题。一阶微分方程的初始条件，如 二阶微分方程的初始条件，如 给出方程解在定义区间两端点处的某种值，称为边界条件，方程和边界条件构成边值问题。00|yyxx,

5、|0000yyyyxxxx及bxayyxfy ),(;)(,)(byay结合下述三种边界条件之一：;)(,)(byay;)()(,)()(1010bybyayay(3)式中。0, 0, 00000它们分别为第一、第二、第三边界问题。方程：(1)(2)(3)n阶初值问题的转化高阶微分方程高阶微分方程( (或方程组或方程组) )的初值问题的初值问题, ,原则上都可以归结为一阶原则上都可以归结为一阶方程组来求解。例如方程组来求解。例如, ,有二阶微分方程的初值问题有二阶微分方程的初值问题 0000)(,)(),(yxyyxyyyxfy在引入新的变量在引入新的变量 后后, ,即化为一阶方程组初值问题即

6、化为一阶方程组初值问题: :yz0000)(,)(,),(yxzyxyzyzyxfz微分方程的解析解与数值解0(0)21yxyyxtxtxy0dee)(22解：但在许多工程实践和科学研究中，除了少数特殊类型的微分方程能用解析法求得其精确解外，大多数情况要得出解的解析表达式是极其困难的，甚至是不可能的。常微分方程的数值问题：由于高阶常微分方程可以转化为一阶常微分方程组，因此，为了不失一般性，本章主要介绍一类一阶常微分方程初值问题：00)()(yxyyx,fy10.110.1 欧拉（欧拉（EulerEuler）法法 10.110.1.1 .1 欧拉法的基本思想欧拉法的基本思想 虽然欧拉法在实际应用

7、时较少采用，但由于它的结构简单，在理论上具有非常重要的意义。欧拉法的基本思想就是用差分方程初值问题 (10.3) 的解来近似微分方程初值问题(10.2)的解，其中 ，式(10.3)也称为欧拉公式。 )()0,1,()(01ayyiy,xfhyyiiii2)-(/abh 欧拉法的几何意义是用一条自点 出发的折线去逼近积分曲线 ，如图10.1所示。因此，这种方法又称为折线法。显然，欧拉法简单地取折线的端点作为数值解，精度非常差。 )(00y,x)(xfy 10.110.1.2 .2 实现欧拉法的基本步骤实现欧拉法的基本步骤（1 1）输入）输入 的区间的区间 ，区间的等分数，区间的等分数 以及以及

8、在在 处处 的初始值的初始值 ；（2 2）对）对 ，计算：，计算： （3 3）输出）输出 ；（4 4）结束。）结束。xnx,x0Ny0 x0y10,1,2,Ni)(iiiy,xff iiifhyyhi*xxiiN/xxhn)-(0Ny例例10.10.1 1 用欧拉法求解初值问题用欧拉法求解初值问题 的数值解（保留的数值解（保留4 4位小数）。位小数）。02(0)y10.y,xxyy# #include include float float funcfunc(float (float x,floatx,float y) y) return(y-x); return(y-x); float fl

9、oat eulereuler(float x0,float (float x0,float xn,floatxn,float y0,int N) y0,int N) float float x,y,hx,y,h; ; intint i i; ; x=x0; x=x0; y=y0; y=y0; h=(xn-x0)/(float)N; / h=(xn-x0)/(float)N; /* * 计算步长计算步长 * */ / for(for(i i=1;i=1;i=N;iN;i+) /+) /* * 欧拉公式欧拉公式 * */ / y=y=y+hy+h* *funcfunc( (x,yx,y);); x

10、=x0+i x=x0+i* *h;h; return(y); return(y); main()main() float x0,xn,y0,e; float x0,xn,y0,e; intint N; N; clrscrclrscr();(); printfprintf(“(“ninputninput n:n ”); / n:n ”); / 输入区间等分数输入区间等分数 scanfscanf(%(%d,&Nd,&N);); printfprintf(“input x0,xn:n ”); / (“input x0,xn:n ”); / 输入输入x x的区间的区间 scanfsca

11、nf(%f%f,&x0,&xn);(%f%f,&x0,&xn); printfprintf(input y0:n ); / (input y0:n ); / 输入输入x x0 0处的处的y y的值的值 scanfscanf(%f,&y0);(%f,&y0); e= e=eulereuler(x0,xn,y0,N);(x0,xn,y0,N); printfprintf(y(%f)=%6.4f,y0,e);(y(%f)=%6.4f,y0,e); 程序运行结果：程序运行结果：input n:input n:5 5input x0,xn:input x0

12、,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.4883y(2.000000)=4.4883input n:input n:1010input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.5937y(2.000000)=4.5937input n:input n:2020input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000

13、)=4.6533y(2.000000)=4.65331e)(xxyx此问题的精确解为此问题的精确解为 ，相应计算值为相应计算值为4.71834.7183。从计算结果。从计算结果可以看出，虽然等分数越大，数值解可以看出，虽然等分数越大，数值解越精确，但精度并不理想。越精确，但精度并不理想。10.1.10.1.3 3 改进的欧拉法改进的欧拉法 10.1.10.1.3 3 改进的欧拉法的基本思想改进的欧拉法的基本思想 改进的欧拉法以隐函数的方式给出 ，计算时常使用迭代法，一般可先由欧拉公式计算出 的初始值 ，再按下面的格式进行 的迭代计算 (10.4) 改进的欧拉法由于添加了迭代过程，计算量较欧拉法

14、增加了一倍，付出这种代价的目的是为了提高精度。1iy1iy(0)1iy1iy)0,1,2,()()(2)()(111)(1(0)1ky,xfy,xfhyyy,xfhyykiiiiikiiiii10.1.410.1.4 实现改进的欧拉法的基本步骤实现改进的欧拉法的基本步骤（1 1）输入）输入 的区间的区间 ，区间的等分数，区间的等分数 以及以及 在在 处处 的初始值的初始值 ；（2 2） , , ， ；（3 3） ；（4 4）计算：）计算：（5 5）若）若 ， ，转（，转（4 4）；否则，转（）；否则，转（6 6）；）；（6 6）输出）输出 ；（7 7）结束。）结束。xnx,x0Ny0 x0y/

15、N)x-x(h0n0 xx 0yy 1i)yx,(fhyyphi*xx0)(pcyx,fh*yy2(/ )yyycpNi 1 iiy# #include include float float funcfunc(float (float x,floatx,float y) y) return(y-x); return(y-x); float float eulereuler(float x0,float (float x0,float xn,floatxn,float y0,int N) y0,int N) float float x,y,yp,yc,hx,y,yp,yc,h; ; intin

16、t i i; ; x=x0; x=x0; y=y0; y=y0; h=(xn-x0)/(float)N; / h=(xn-x0)/(float)N; /* * 计算步长计算步长 * */ / for(for(i i=1;i=1;i=N;iN;i+) /+) /* * 改进的欧拉公式改进的欧拉公式 * */ / ypyp= =y+hy+h* *funcfunc( (x,yx,y);); x=x0+i x=x0+i* *h;h; ycyc= =y+hy+h* *funcfunc( (x,ypx,yp);); y=( y=(yp+ycyp+yc)/2.0;)/2.0; return(y); retu

17、rn(y); 例例10.10.1 1 用改进的欧拉法求解初值问题用改进的欧拉法求解初值问题 的数值解（保留的数值解（保留4 4位小数）。位小数）。02(0)y10.y,xxyy程序运行结果：程序运行结果：input n:input n:5 5input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.7027y(2.000000)=4.7027input n:input n:1010input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0

18、: 2.0 2.0y(2.000000)=4.7141y(2.000000)=4.7141input n:input n:2020input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.7172y(2.000000)=4.717210.1.5 欧拉方法的收敛性 局部截断误差由欧拉公式计算一步得到的误差记为： 11212 () Taylor, (,),()()()()()2 ()(, ()() (1)2nnnnnnnnnnnnnny xxxxhy xy xhy xhy xyhy xhf xy

19、xy将在点展开1(),(8.2) () (, () (2)nnnnny xyy xh f xy x用准确值 利用欧拉公式计算得到的称为“局部截断误差”。令 定义 若给定方法的局部截断误差满足则称该方法是 p 阶的，或称为具有 p 阶精度。11112() ()()(, () (). 2nnnnnnnnTy xyy xy xhf xy xhy2221 max |( )|,|()|(). 22a x bnnMy xhhTyMO h 则11|(),pnTO h 整体截断误差 记111().nnney xy112121112 ,(), (), () ,(),.nnnnnnyy yyy xy xy xy

20、xey yy因为计算 时 用到的, , 是 的近似值 每步产生的误差会累积到计算的误差中因此 与 , , 都有关 称为整体截断误差1111111111| | ()| | ()| |. (3)nnnnnnnnnney xyy xyyyTyy1111,(8.2)| | () (, () (,) | | () |(, ()(,)| | () | () | () (1)nnnnnnnnnnnnnnnnnnnnyyyyy xh f xy xyh f xyy xyh f xy xf xyy xyhL y xyhL对应用欧拉公式得李普希兹条件|. (4)ne由此知，当 max |,(4)(3)kkTT记 将

21、代入得1121210112|(1)| (1)(1)| (1)(1) | (1)(1)(1) (1)|(1)1(1)1 ()(1) 1 (nnnnnnnneThLeThL ThLeThL ThLeThL ThL ThL ThLehLhLTO hhLhLO).h10,|0, nhe有欧拉公式是一阶收敛的。 注 整体截断误差与局部截断误差的关系： 类似于向前欧拉公式收敛性的分析，可以得到： 向后欧拉公式： 局部截断误差： 整体截断误差： 收敛的条件： 1(1)() 0,(1)e , (1)ee.nnznnhLL b azzhLconst111|(|).nneO hT21|(),nTO h1|( ).

22、neO h01,8.4hL因为由（）知，.(1)( )( )(1)( )(1)11111111|(,)(,)|.Lkkkkkknnnnnnnnyyh f xyf xyhL yy 梯形公式：（二阶方法） 局部截断误差： 整体截断误差： 收敛的条件： 31|();nTO h21|().neO h101,8.92hL因为由（） 知，(1)( )( )(1)111111.( )(1)11|(,)(,) |.2kkkknnnnnnLkknnyyh f xyf xyhLyy10.2 龙格龙格- -库塔库塔(Runge-Kutta)方法方法 回顾： 向前Euler: 向后Euler: 梯形公式: 研究提高收

23、敛阶的方法 用Taylor展式|( );neO h|( );neO h2|().neO h21( )1()()()()()2! ()(),!nnnnnpppnhy xy xhy xhy xy xhyxO hp采用此式，需要计算高阶导数，如112( )1 (), .2!nnppnnnnny xyhhyyhyyyp22( , ),( , ),()2, xyxyxyxxxyxyyyyyf x yyfx yffyfffyfffffff ffff f 10.2.1 二阶二阶Runge-Kutta)方法方法 基本思想：分别计算不同点上的函数值，然后对这些函数值进行线性组合，构造出近似计算公式，然后把此近似

24、公式与解 的Taylor展开式进行比较，使他们前面的若干项保持一致，从而使近似公式达到需要的收敛阶。 考察梯形公式： 引入参数 k1， k2 ，上式改写为1()ny x111 (,)(,) .2nnnnnnhyyf xyf xy 为便于后面的推广，将上式写成更一般的形式 任务：适当选取参数 使公式（2）的收敛阶尽可能的高。11212111 2(,) (1) (,)(,) nnnnnnnnhyykkkf xykf xyf xh yhk11 122122211 (,) (2) (,) nnnnnnyyh c kc kkf xykf xa h yb hk12221,c c a b 先将k2在(xn,

25、 yn)处Taylor展开 将k1 ，k2 代入（2)中的第一式，得23123()()()()()2! ()(, ()(, ()2 (, ()(, ()() (4)nnnnnnnxnnynnnnhy xy xhy xy xO hhy xhf xy xfxy xfxy xf xy xO h而 22211222112221(,) (,)(,)(,)() (,)(,)(,)(,)()nnnnxnnynnnnxnnnnynnkf xa h yb hkf xya hfxyb hk fxyO hf xya hfxyb hf xyfxyO h2112222321 2(,)(,) (,)(,)() (3)nn

26、nnxnnnnynnyyh ccf xya c h fxyb c h f xyfxyO h要使(3)逼近(4)具有O(h3) 的截断误差,比较两式的系数得 特别地， 取12222211 2121cca ca b122211, 1,(2)2ccab 由得梯形公式：112121() 2(,) (8.15) (,) nnnnnnhyykkkf xykf xh yhk 取12121 (,) (8.16) (,)22nnnnnnyyhkkf xyhhkf xyk12221101, ,(2)2ccab，由即得中点公式：(5),(6)-也称为二阶龙格 库塔公式.2（计算 个点函数值）3121|(), |()

27、.nnTO heO h 10.2.2 三阶三阶Runge-Kutta方法方法 增加一个参与运算的函数值，设 将k2 , k3在(xn, yn)处Taylor展开到三阶，然后将k2 , k2 , k3代入（5）式，合并同类项，然后与 y(xn+1) 在 xn 处Taylor展开对应项相比较，得参数满足的方程：11 1223312221133311322 (, ) (5)(, ) (, ) nnnnnnnnyyh c kc kc kkf xykf xa h yb hkkf xa h yb hkb hk12322133132223322223332 321 1 21 31 6cccababbc ac

28、 ac ac ac a b (6) 取满足（6）式，代入（5）式，得三阶R-K公式（10.17）：12322133132121, ,6361, 1, 1, 2,2cccababb ，11231213124 6(, ) (8.17)11(, )22 (, 2) nnnnnnnnhyykkkkf xykf xh yhkkf xh yhkhk 取满足（6）式，代入（5）式，得三阶R-K公式（10.18）：12322133132130, ,44122, , 0, ,333cccababb，113121323 4(, ) (8.18)11(, )3322 (, ) 33nnnnnnnnhyykkkf x

29、ykf xh yhkkf xh yhk4131|()|()nnTO heO h 取满足（6）式，代入（5）式，得三阶R-K公式（10.19）：12322133132214, ,939133, , 0, ,244cccababb，112312132234 9(, ) (8.19)11(, )2233 (, ) 44nnnnnnnnhyykkkkf xykf xh yhkkf xh yhk10.2.3 四阶四阶Runge-Kutta方法方法11234121324322 6(, ) 11(, ) 2211(, ) 22 (, ) nnnnnnnnnnhyykkkkkf xykf xh yhkkf x

30、h yhkkf xh yhk(8.20) 局部截断误和整体截断误差：11234121312412333 8(, ) 11(, ) 3321(, ) 33(, ) nnnnnnnnnnhyykkkkkf xykf xh yhkkf xh yhkhkkf xh yhkhkhk(8.21)5141|()|()nnTO heO h10.10.4.2 4.2 实现标准四阶龙格实现标准四阶龙格库塔法的基本步骤库塔法的基本步骤（1 1）输入）输入 的区间的区间 ，区间的等分数，区间的等分数 以及以及 在在 处的初始值处的初始值 ；（2 2） , , ， ；（3 3） ；（4 4）计算：）计算：（5 5）若）

31、若 ， ，转（，转（4 4）；否则，转（）；否则，转（6 6）；）；（6 6）输出）输出 ；（7 7）结束。）结束。xnx,x0Ny0 x0y/Nxxhn)-(00 xx 0yy 1i2/hxxh)(1yx,fd )21(2/dh*y,xfdh)22(3/dh*y,xfdh)3(4dh*y,xfdh4)/63*22*21(ddddh*yyhi*xx0Ni 1 iiy例例10.10.3 3 用标准四阶龙格用标准四阶龙格库塔法求解初值问题的数值库塔法求解初值问题的数值 解（保留解（保留6 6位小数）。位小数）。0(0)10yxyxy#include #include float float fun

32、cfunc(float (float x,floatx,float y) y) return(x-y); return(x-y); float float runge_kuttarunge_kutta(float x0,float (float x0,float xn,floatxn,float y0,int N) y0,int N) float x,y,y1,y2,h,xh; float x,y,y1,y2,h,xh; float d1,d2,d3,d4; float d1,d2,d3,d4; intint i i; ; x=x0; x=x0; y=y0; y=y0; h=(xn-x0)/(

33、float)N; / h=(xn-x0)/(float)N; /* * 计算步长计算步长 * */ / for(for(i i=1;i=1;i=N;iN;i+) /+) /* * 标准四阶龙格标准四阶龙格- -库塔公式库塔公式 * */ / xhxh= =x+hx+h/2;/2; d1= d1=funcfunc( (x,yx,y);); d2= d2=funcfunc( (xh,y+hxh,y+h* *d1/2.0);d1/2.0); d3= d3=funcfunc( (xh,y+hxh,y+h* *d2/2.0);d2/2.0); d4= d4=funcfunc( (xh,y+hxh,y+h* *d3);d3); y= y=y+hy+h* *(d1+2(d1+2* *d2+2d2+2* *d3+d4)/6.0;d3+d4)/6.0; x=x0+i x=x0+i* *h; h; return(y); return(y); main()main() float x0,xn,y0,e; float x0,xn,y0,e; intint N; N; clrscrclrscr();(); printfprintf(ninputninput n:n ); / n:n ); /输入区间等分数输入区间等分数 scanfscanf(%(%