问题描述
我正在编写代码来使用 Dormand Prince 方法求解微分方程组。我在堆栈溢出中找到了相同的 C 代码，并尝试为我的方程实现它。它起作用了，我得到了解决方案。但是，我需要通过 fortran 和 C 来解决它。当我将代码转换为 fortran 时，它只会给出无限次相同的结果。我是 Fortran 的初学者。请帮助我哪里出错了？`
#include<stdio.h>
#include<math.h>
#include<stdbool.h>
#include <stdlib.h>
#define N 3
double ddopri5(void fcn(double,double *,double *),double *y);
double alpha;
void fcn(double t,double *y,double *f);
double eps;
int main(void){
double y[N];
//eps = 1.e-9;
printf("Enter alpha:\n");
scanf("%lg",&alpha);
printf("Enter epsilon:\n");
scanf("%lg",&eps);
y[0]=1.0;//x1(0)
y[1]=1.0;//x2(0)
y[2]=1.0;//p1(0)
ddopri5(fcn,y);
}
void fcn(double t,double *f){
/* double h = 0.25;*/
f[0] = y[1];
f[1] =- y[2]*y[0];
f[2] = t;
}
double ddopri5(void fcn(double,double *y){
double t,h,a,b,tw,chi;
double w[N],k1[N],k2[N],k3[N],k4[N],k5[N],k6[N],k7[N],err[N],dy[N];
int i;
double errabs;
int iteration;
iteration = 0;
//eps = 1.e-9;
h = 0.5;
a = 0.0;
b = 15.0;//3.1415926535;
t = a;
while(t < b -eps){
/* printf("%lg\n",eps);*/
fcn(t,y,k1);
tw = t+ (1.0/5.0)*h;
for(i = 0; i < N; i++){
/*printf("k1[%i] = %.15lf \n",i,k1[i]);*/
w[i] = y[i] + h*(1.0/5.0)*k1[i];
}
fcn(tw,w,k2);
tw = t+ (3.0/10.0)*h;
for(i = 0; i < N; i++){
/*printf("k2[%i] = %.15lf \n",k2[i]);*/
w[i] = y[i] + h*((3.0/40.0)*k1[i] + (9.0/40.0)*k2[i]);
}
fcn(tw,k3);
tw = t+ (4.0/5.0)*h;
for(i = 0; i < N; i++){
/*printf("k3[%i] = %.15lf \n",k3[i]);*/
w[i] = y[i] + h*((44.0/45.0)*k1[i] - (56.0/15.0)*k2[i] + (32.0/9.0)*k3[i]);
}
fcn(tw,k4);
tw = t+ (8.0/9.0)*h;
for(i = 0; i < N; i++){
/*printf("k4[%i] = %.15lf \n",k4[i]);*/
w[i] = y[i] + h*((19372.0/6561.0)*k1[i] - (25360.0/2187.0)*k2[i] +
(64448.0/6561.0)*k3[i] - (212.0/729.0)*k4[i]);
}
fcn(tw,k5);
tw = t + h;
for(i = 0; i < N; i++){
/*printf("k5[%i] = %.15lf \n",k5[i]);*/
w[i] = y[i] + h*((9017.0/3168.0)*k1[i] - (355.0/33.0)*k2[i] +
(46732.0/5247.0)*k3[i] + (49.0/176.0)*k4[i] - (5103.0/18656.0)*k5[i]) ;
}
fcn(tw,k6);
tw = t + h;
for(i = 0; i < N; i++){
/*printf("k6[%i] = %.15lf \n",k6[i]);*/
w[i] = y[i] + h*((35.0/384.0)*k1[i] + (500.0/1113.0)*k3[i] +
(125.0/192.0)*k4[i] - (2187.0/6784.0)*k5[i] + (11.0/84.0)*k6[i]);
}
fcn(tw,k7);
errabs = 0;
for(i = 0; i < N; i++){
/* printf("k7[%i] = %.15lf \n",k7[i]);*/
/* dy[i] = h*((71.0/57600.0)*k1[i] - (71.0/16695.0)*k3[i] +
(71.0/1920.0)*k4[i] - (17253.0/339200.0)*k5[i] + (22.0/525.0)*k6[i]);*/
dy[i] = h*((35.0/384.0)*k1[i] + (500.0/1113.0)*k3[i] + (125.0/192.0)*k4[i] -
(2187.0/6784.0)*k5[i] + (11.0/84.0)*k6[i]);
/*err[i] = h*((71.0/57600.0)*k1[i] + (71.0/16695.0)*k3[i] +
(71.0/1920.0)*k4[i] - (17253.0/339200.0)*k5[i] + (22.0/525.0)*k6[i] -
(1.0/40.0)*k7[i])*/;
err[i] = h*((71.0/57600.0)*k1[i] - (71.0/16695.0)*k3[i] +
(71.0/1920.0)*k4[i] - (17253.0/339200.0)*k5[i] + (22.0/525.0)*k6[i] -
(1.0/40.0)*k7[i]);
/*printf("err[%i] = %.15lf \n",err[i]);*/
errabs+=err[i]*err[i];
}
errabs = sqrt(errabs);
/* printf("errabs = %.15lf\n",errabs); */
if( errabs < eps){
t+= h;
/* printf(" FROM IF \t t = %.25lf,\n h = %.25lf,\n errabs = %.25lf,\n
iteration = %i . \n",t,errabs,iteration);*/
for(i = 0; i < N; i++){
y[i]+=dy[i];
}
}
chi=errabs/eps;
chi = pow(chi,(1.0/6.0));
if(chi > 10) chi = 10;
if(chi < 0.1) chi = 0.1;
h*= 0.95/chi;
if( t + h > b ) h = b - t;
/* for(i = 0; i < N; i++){
printf("y[%i] = %.15lf \n",y[i]);
}*/
iteration++;
/*printf("t = %.25lf \t h = %.25lf\n",h);*/
/*if(iteration > 5) break;*/
/* printf("end \n");*/
for(i = 0; i < N; i++){
printf("%.15lf %.15lf %.15lf %.15lf\n",y[0],y[1],y[2]);
}
if(iteration > 30000) break;
}
/* for(i = 0; i < N; i++){
printf("y[%i] = %.15lf\n",y[i]);
}*/
return 0;
}
转换后的Fortran版本如下
PROGRAM GHO
IMPLICIT NONE
DOUBLE PRECISION::y(3),eps,f(3)
y(1)=1.0
y(2)=-1.22565282791
y(3)=-0.274772807644
eps=0.05
call ddopri5( y)
END PROGRAM
SUbroUTINE fcn(t,f)
IMPLICIT NONE
DOUBLE PRECISION::y(3),f(3)
f(1)=y(2)
f(2)=-y(3)*y(1)
f(3)=t
END SUbroUTINE
SUbroUTINE ddopri5( y)
IMPLICIT NONE
external fcn
DOUBLE PRECISION::t,chi,eps
DOUBLE PRECISION:: w(3),k1(3),k2(3),k3(3),k4(3),k5(3),k6(3),k7(3),&
yerr(3),dy(3),y(3)
INTEGER :: i
DOUBLE PRECISION::errabs
INTEGER::iteration;
iteration = 0
h=0.5
a=0.0
b=15.0
t=a
do while(t.le.(b-eps))
call fcn(t,k1)
tw = t+ (1.0/5.0)*h
print *,k1(1),k1(2),k1(3)
do i=1,3
w(i) = y(i)+ h*(1.0/5.0)*k1(i)
enddo
call fcn(tw,k2)
tw = t+ (3.0/10.0)*h
do i=1,3
w(i) = y(i)+ h*((3.0/40.0)*k1(i)+ (9.0/40.0)*k2(i))
enddo
call fcn(tw,k3)
tw = t+ (4.0/5.0)*h
do i=1,3
w(i) = (i)+ h*((44.0/45.0)*k1(i)- (56.0/15.0)*k2(i)+ (32.0/9.0)*k3(i))
enddo
call fcn(tw,k4)
tw = t+ (8.0/9.0)*h
do i=1,3
w(i) = y(i)+ h*((19372.0/6561.0)*k1(i)- (25360.0/2187.0)*k2(i) &
+(64448.0/6561.0)*k3(i)- (212.0/729.0)*k4(i))
enddo
call fcn(tw,k5)
tw = t + h;
do i=1,3
w(i)= y(i)+h*((9017.0/3168.0)*k1(i)-(355.0/33.0)*k2(i)&
+(46732.0/5247.0)*k3(i)&
+ (49.0/176.0)*k4(i)-(5103.0/18656.0)*k5(i))
enddo
call fcn(tw,k6)
tw = t + h
do i=1,3
w(i)= y(i)+ h*((35.0/384.0)*k1(i)+(500.0/1113.0)*k3(i)+
(125.0/192.0)*k4(i)&
-(2187.0/6784.0)*k5(i)+(11.0/84.0)*k6(i))
enddo
call fcn(tw,k7)
errabs = 0
do i=1,3
dy(i)= h*((35.0/384.0)*k1(i)+ (500.0/1113.0)*k3(i)+ (125.0/192.0)*k4(i)&
- (2187.0/6784.0)*k5(i)+ (11.0/84.0)*k6(i))
yerr(i)=h*((71.0/57600.0)*k1(i)-(71.0/16695.0)*k3(i)+(71.0/1920.0)*k4(i)&
-(17253.0/339200.0)*k5(i)+(22.0/525.0)*k6(i)-(1.0/40.0)*k7(i))
print *,yerr(1),yerr(2),yerr(3)
errabs=errabs+MAXVAL((yerr(i)*yerr(i)))
print *,errabs
enddo
errabs =sqrt(errabs)
print *,errabs
if( errabs<eps)then
t=t+h
do i=1,3
y(i)=y(i)+dy(i)
enddo
endif
chi=errabs/eps
chi =chi**(1.0/6.0)
if(chi>10) chi = 10
if(chi<0.1) chi = 0.1
h=(0.95/chi)*h
print *,h
if( t + h > b ) h=b-t
iteration=iteration+1
write (*,*) t,y(1),y(2),y(3)
enddo
END SUbroUTINE
请帮帮我。我被困在这里好几天了。
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）