24小时热门版块排行榜

返回列表

yuwuxiaoyuan

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 444.8
帖子: 50
在线: 13.4小时
虫号: 2382452
注册: 2013-03-27
专业: 工程热物理相关交叉领域

[求助] c语言龙阁库塔问题求助各位大神

在建模计算的过程中用到了四阶龙阁库塔方法，然后我改编了程序集中的一段程序，可是各种问题。要么就是赋完初值之后，后面的计算结果一直都是初值。因为老师要求的特别急，所以特此求助各位大神～
这个是main函数：
#include "stdio.h"
#include "grkt2.c"
  main()
  { int i;
double t,h,eps,y[4];
void grkt2f();
y[0]=0.0075;y[1]=0.1;y[2]=25;y[3]=23;
t=0.0; h=0.1; eps=0.0001;
printf("\n" );
printf("t=%7.3f  y(0)=%e  y(1)=%e  y(2)=%e  y(3)=%e\n",t,y[0],y[1],y[2],y[3]);
for (i=1; i<=500; i++)
   { grkt2(t,h,y,4,eps,grkt2f);
      t=t+h;
      printf("t=%7.3f  y(0)=%e  y(1)=%e  y(2)=%e  y(3)=%e\n",t,y[0],y[1],y[2],y[3]);
   }
  }

  void grkt2f()
{
int n;
double t,y[4],d[4],A[4],B[4],C3,C4,k[4],a,b,e,f,s,v,aln,alq,Aln,Aghw,Tln,Gt,Ghw,GH,Mg,Mhw,Cg,Chw,Thwc,Thwr;
t=t; n=n;
s=100.406;    /*冷凝器底面积*/
v=1348.9546;    /*冷凝器箱体体积*/
aln=6.3735;    /*冷凝水液膜换热系数*/
alq=29.17134; /*循环水与管壁换热系数*/
Aln=17883;    /*冷凝管与蒸汽的换热面积*/
Aghw=17800;    /*冷凝管与循环水的换热面积*/
Tln=40.23;    /*冷凝器内平均温度（摄氏度）*/
Gt=276.742;    /*冷凝器内净增加的质量流量*/
Ghw=44960;    /*循环水质量流量*/
GH=2041574.281;  /*冷凝器内净增加的焓值*/
Mg=2816572.5; /*钛合金管材总质量*/
Mhw=113551.144;  /*管内循环水的总质量*/
Cg=0.52/1000; /*钛合金管材的比热容*/
Chw=4.174;    /*循环水的定压比热容*/
Thwc=33.300;    /*循环水出口温度*/
Thwr=23.000;    /*循环水进口温度*/
a=22833636.71;b=1098193947;e=8744.2610;f=-650.8568;
k[1]=999.95;k[2]=e+f*y[2];k[3]=168516.21;k[4]=a+b*y[2];
A[1]=(Gt*k[4]-GH*k[2]+Aln*Tln*aln*k[2])/(s*(k[1]*k[4]-k[2]*k[3]));
A[2]=(Gt*k[3]-GH*k[1]+Aln*Tln*aln*k[1])/(k[2]*k[3]-k[1]*k[4]);
A[3]=Aln*Tln*aln/(Mg*Cg);
A[4]=-Ghw*Chw*(Thwc-Thwr)/(Mhw*Chw);
B[1]=-k[2]*Aln*aln/(s*(k[1]*k[4]-k[2]*k[3]));
B[2]=-k[1]*Aln*aln/(k[2]*k[3]-k[1]*k[4]);
B[3]=-(aln*Aln+alq*Aghw)/(Mg*Cg);
B[4]=alq*Aghw/(Mhw*Chw);
C3=Aghw*alq/(Mg*Cg);
C4=-alq*Aghw/(Mhw*Chw);
d[0]=A[1]+B[1]*y[2];
d[1]=A[2]+B[2]*y[2];
d[2]=A[3]+B[3]*y[2]+C3*y[3];
d[3]=A[4]+B[4]*y[2]+C4*y[3];
return;
}

然后是调用的龙阁库塔的函数，这段是直接在程序集上的（徐士良版）

#include "stdlib.h"
  #include "math.h"
  void grkt2(t,h,y,n,eps,grkt2f)
void grkt2f()

int n;
double t,h,eps,y[];
{int m,i,j,k;
double hh,p,dt,x,tt,q,a[4],*g,*b,*c,*d,*e;
g=malloc(n*sizeof(double));
b=malloc(n*sizeof(double));
c=malloc(n*sizeof(double));
d=malloc(n*sizeof(double));
e=malloc(n*sizeof(double));
hh=h; m=1; p=1.0+eps; x=t;
for (i=0; i<=n-1; i++) c=y;
while (p>=eps)
   { a[0]=hh/2.0; a[1]=a[0]; a[2]=hh; a[3]=hh;
      for (i=0; i<=n-1; i++)
      { g=y; y=c;}
      dt=h/m; t=x;
      for (j=0; j<=m-1; j++)
      { grkt2f(t,y,n,d);
         for (i=0; i<=n-1; i++)
            { b=y; e=y;}
         for (k=0; k<=2; k++)
            { for (i=0; i<=n-1; i++)
               { y=e+a[k]*d;
                  b=b+a[k+1]*d/3.0;
               }
            tt=t+a[k];
            grkt2f(tt,y,n,d);
            }
         for (i=0; i<=n-1; i++)
            y=b+hh*d/6.0;
         t=t+dt;
      }
      p=0.0;
      for (i=0; i<=n-1; i++
      { q=fabs(y-g);
         if (q>p) p=q;
      }
      hh=hh/2.0; m=m+m;
   }
free(g); free(b); free(c); free(d); free(e);
return;
  }

[ Last edited by yuwuxiaoyuan on 2013-4-8 at 12:53 ]