常微分方程数值解rk4ab4am4.cpp
来自「数值计算各算法的程序,包括:牛顿迭代法,超松驰迭代法,微分法等.对于初学数值分析」· C++ 代码 · 共 234 行
CPP
234 行
#include <iostream.h>
#include <math.h>
#include <iomanip.h>
double f(double x,double y)
{
double f;
f=-pow(x,2)*pow(y,2);
return f;
}
double exactvalue(double x)
{
double e;
e=3/(1+pow(x,3));
return e;
}
void RK4(double y0,double h,double a,int n,double *x)
{
double e,err,k1,k2,k3,k4;
int i;
double *y;
y=new double[n];
y[0]=e=y0;
for(i=0;i<n-1;i++)
{
k1=f(x[i],y[i]);
k2=f(x[i]+h/2,y[i]+h*k1/2);
k3=f(x[i]+h/2,y[i]+h*k2/2);
k4=f(x[i]+h,y[i]+h*k3);
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
}
cout<<"用RK4方法计算："<<endl;
cout<<"计算结果 精确解 误差"<<endl;
cout<<y[0]<<setw(20)<<e<<setw(20)<<fabs(y[0]-e)<<endl;
for(i=0;i<n-1;i++)
{
e=exactvalue(x[i+1]);
err=fabs(y[i+1]-e);
cout<<y[i+1]<<setw(20)<<e<<setw(20)<<err<<endl;
}
delete y;
}
void AB4(double y0,double h,int n,double *x)
{
int i;
double e=y0,err,k1,k2,k3,k4;
double *y;
y=new double[n];
y[0]=y0;
for(i=0;i<=2;i++)
{
k1=f(x[i],y[i]);
k2=f(x[i]+h/2,y[i]+h*k1/2);
k3=f(x[i]+h/2,y[i]+h*k2/2);
k4=f(x[i]+h,y[i]+h*k3);
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
}
for(i=3;i<n-1;i++)
{
y[i+1]=y[i]+h*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])
-9*f(x[i-3],y[i-3]))/24;
}
cout<<"用AK4方法计算："<<endl;
cout<<"计算结果 精确解 误差"<<endl;
cout<<y[0]<<setw(20)<<e<<setw(20)<<fabs(y[0]-e)<<endl;
for(i=0;i<n-1;i++)
{
e=exactvalue(x[i+1]);
err=fabs(y[i+1]-e);
cout<<y[i+1]<<setw(20)<<e<<setw(20)<<err<<endl;
}
}
void AB4AM4(double y0,double h,int n,double *x)
{
int i;
double k1,k2,k3,k4,e=y0,err;
double *y;
y=new double[n];
y[0]=y0;
for(i=0;i<=2;i++)
{
k1=f(x[i],y[i]);
k2=f(x[i]+h/2,y[i]+h*k1/2);
k3=f(x[i]+h/2,y[i]+h*k2/2);
k4=f(x[i]+h,y[i]+h*k3);
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
}
for(i=3;i<n-1;i++)
{
y[i+1]=y[i]+h*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])
-9*f(x[i-3],y[i-3]))/24;
y[i+1]=y[i]+h*(9*f(x[i+1],y[i+1])+19*f(x[i],y[i])
-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]))/24;
}
cout<<"用AK4AM4预测校正方法计算："<<endl;
cout<<"计算结果 精确解 误差"<<endl;
cout<<y[0]<<setw(20)<<e<<setw(20)<<fabs(y[0]-e)<<endl;
for(i=0;i<n-1;i++)
{
e=exactvalue(x[i+1]);
err=fabs(y[i+1]-e);
cout<<y[i+1]<<setw(20)<<e<<setw(20)<<err<<endl;
}
}
void AB4AM4p(double y0,double h,int n,double *x)
{
int i;
double k1,k2,k3,k4,e=y0,err,temp;
double *y;
y=new double[n];
y[0]=y0;
for(i=0;i<=2;i++)
{
k1=f(x[i],y[i]);
k2=f(x[i]+h/2,y[i]+h*k1/2);
k3=f(x[i]+h/2,y[i]+h*k2/2);
k4=f(x[i]+h,y[i]+h*k3);
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
}
for(i=3;i<n-1;i++)
{
y[i+1]=y[i]+h*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])
-9*f(x[i-3],y[i-3]))/24;
temp=y[i+1];
y[i+1]=y[i]+h*(9*f(x[i+1],y[i+1])+19*f(x[i],y[i])
-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]))/24;
y[i+1]=251*y[i+1]/270+19*temp/270;
}
cout<<"用AK4AM4预测校正方法计算："<<endl;
cout<<"计算结果 精确解 误差"<<endl;
cout<<y[0]<<setw(20)<<e<<setw(20)<<fabs(y[0]-e)<<endl;
for(i=0;i<n-1;i++)
{
e=exactvalue(x[i+1]);
err=fabs(y[i+1]-e);
cout<<y[i+1]<<setw(20)<<e<<setw(20)<<err<<endl;
}
}
void main()
{
double a=0,b=1.5,h=0.1,y0=3;
int n;
n=int((b-a)/h)+1;
int i;
double *x;
x=new double[n];
for(i=0;i<=n-1;i++)
{
x[i]=a+h*i;
}
RK4(y0,h,a,n,x);
AB4(y0,h,n,x);
AB4AM4(y0,h,n,x);
AB4AM4p(y0,h,n,x);
}
/*程序运行结果为：
用RK4方法计算：
计算结果 精确解 误差
3 3 0
2.997 2.997 1.87138e-007
2.97619 2.97619 3.91665e-007
2.92113 2.92113 7.58342e-007
2.81955 2.81955 1.61101e-006
2.66666 2.66667 3.17735e-006
2.4671 2.46711 5.00551e-006
2.2338 2.2338 5.77233e-006
1.98412 1.98413 4.12954e-006
1.73511 1.73511 1.15554e-007
1.50001 1.5 5.80668e-006
1.28701 1.287 1.13075e-005
1.09972 1.09971 1.54242e-005
0.938397 0.93838 1.77272e-005
0.8013 0.801282 1.83754e-005
0.685732 0.685714 1.78e-005
用AK4方法计算：
计算结果 精确解 误差
3 3 0
2.997 2.997 1.87138e-007
2.97619 2.97619 3.91665e-007
2.92113 2.92113 7.58342e-007
2.81839 2.81955 0.00115961
2.66467 2.66667 0.0019942
2.4652 2.46711 0.00190263
2.23308 2.2338 0.000725963
1.98495 1.98413 0.000823594
1.73704 1.73511 0.00193629
1.50219 1.5 0.00219455
1.28876 1.287 0.00176216
1.10072 1.09971 0.00101746
0.93871 0.93838 0.000330765
0.801135 0.801282 0.000147103
0.685335 0.685714 0.00037971
用AK4AM4预测校正方法计算：
计算结果 精确解 误差
3 3 0
2.997 2.997 1.87138e-007
2.97619 2.97619 3.91665e-007
2.92113 2.92113 7.58342e-007
2.81968 2.81955 0.00012956
2.66688 2.66667 0.000209316
2.46725 2.46711 0.000146501
2.23373 2.2338 7.34996e-005
1.98379 1.98413 0.00034028
1.73461 1.73511 0.000499563
1.49952 1.5 0.000484064
1.28666 1.287 0.000344148
1.09953 1.09971 0.000173599
0.938343 0.93838 3.72081e-005
0.801327 0.801282 4.5323e-005
0.685796 0.685714 8.18254e-005
用AK4AM4预测校正方法计算：
计算结果 精确解 误差
3 3 0
2.997 2.997 1.87138e-007
2.97619 2.97619 3.91665e-007
2.92113 2.92113 7.58342e-007
2.81959 2.81955 3.88409e-005
2.66671 2.66667 4.61875e-005
2.4671 2.46711 8.23448e-006
2.23368 2.2338 0.000122422
1.98388 1.98413 0.0002423
1.73481 1.73511 0.000298987
1.49973 1.5 0.000268093
1.28682 1.287 0.000180611
1.09962 1.09971 8.49646e-005
0.938367 0.93838 1.2415e-005
0.801311 0.801282 2.93031e-005
0.68576 0.685714 4.61604e-005
*/
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?