📄 常微分方程数值解rk4ab4am4.cpp
字号:
#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 + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -