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

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

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

（1-1）

，

（1-2）

（1-3）

（1-4）

（1-5）

（1-6）

（1-7）

（1-8）

（1-9）

（1-10）

经过循环计算由推得……

每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计

算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常

精准，稳定，且易于编程。

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<

cout<< initial[0]<

for(i=0;i

{ 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 经典四阶龙格库塔法解一阶微分方程算法程序调试图2.高斯列主元法解线性方程组

2.1高斯列主元法解线性方程组算法分析

使用伪代码编写高斯消元过程：

for k=1 to n-1 do

for i=k+1 to n

l<=a(i,k)/a(k,k)

for j=k to n do

a(i,j)<=a(i,j)-l*a(k,j)

end %end of for j

b(i)<=b(i)-l*b(k)

end %end of for i

end %end of for k

最后得到A，b可以构成上三角线性方程组

接着使用回代法求解上三角线性方程组

因为高斯消元要求a（k,k）≠0（k=1,2,3……n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为：

①找主元：计算,并记录其所在行r ,

②交换第r行与第k行；

③以第k行为工具行处理以下各行，使得从第k列的第k+1行到第n行的元素全部为0；

④得到增广矩阵的上三角线性方程组；

⑤使用回代法对上三角线性方程组进行求解

2.2高斯列主元法解线性方程组流程图

图2-1 高斯列主元法解线性方程组流程图

2.3高斯列主元法解线性方程组程序代码

#include

#include

#define N 3

using namespace std;

void main()

{int i,j,k,n,p;

float t,s,m,a[N][N],b[N],x[N];

cout<<"请输入方程组的系数"<

for(i=0;i

{for(j=0;j

cin>>a[i][j];}

cout<<"请输入方程组右端的常数项:"<

for(i=0;i

cin>>b[i];

for(j=0;j