经典四阶龙格库塔法解一阶微分方程组

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

1.经典四阶龙格库塔法解一阶微分方程组

1.1运用四阶龙格库塔法解一阶微分方程组算法分析

),,(1k k k y x t f f =, )2,2,2(112g h

y f h x h t f f k k k +++=

)2

,2,2(223g h

y f h x h t f f k k k +++=

),,(334hg y hf x h t f f k k k +++=

),,(1k k k y x t g g =

)2,2,2(112g h

y f h x h t g g k k k +++=

)2,2,2(223g h

y f h x h t g g k k k +++=

),,(334hg y hf x h t g g k k k +++=

)

22(6

)22(6

43211

43211g g g g h

y y f f f f h

x x k k k k ++++=++++=++

1k k t t h +=+

经过循环计算由 推得 ……

每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为()

N O h ,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精

准,稳定,且易于编程。

000,,t x y ()()111222,,,,t x y t x y (1-1)

(1-2)

(1-3)

(1-4)

(1-5)

(1-6)

(1-7)

(1-8)

(1-9)

(1-10)

1.2经典四阶龙格库塔法解一阶微分方程流程图

图1-1 经典四阶龙格库塔法解一阶微分方程流程图

1.3经典四阶龙格库塔法解一阶微分方程程序代码:

#include

#include

using namespace std;

void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h) {

double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;

t0=initial[0];x0=initial[1];y0=initial[2];

f1=f(t0,x0,y0); g1=g(t0,x0,y0);

f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);

f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2,

x0+h*f2/2,y0+h*g2/2);

f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;

}

int main()

{

double f(double t,double x, double y);

double g(double t,double x, double y);

double initial[3],resu[3];

double a,b,H;

double t,step;

int i;

cout<<"输入所求微分方程组的初值t0,x0,y0:";

cin>>initial[0]>>initial[1]>>initial[2];

cout<<"输入所求微分方程组的微分区间[a,b]:";

cin>>a>>b;

cout<<"输入所求微分方程组所分解子区间的个数step:";

cin>>step;

cout<

H=(b-a)/step;

cout<< initial[0]<

{ RK4( f,g ,initial, resu,H);

cout<

initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];

}

return(0);

}

double f(double t,double x, double y)

{

double dx;

dx=x+2*y; return(dx); }

double g(double t,double x, double y) {

double dy; dy=3*x+2*y; return(dy); }

1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示:

应用所编写程序计算所给例题:

其中初值为

求解区间为[0,0.2]。

计算结果为:

图1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图

⎪⎪⎩⎪⎪⎨

⎧+=+=y

x dt

dy y x dt dx

232⎩⎨⎧==4

)0(6)0(y x

相关文档
最新文档