| 查看: 1115 | 回复: 5 | ||
yuwuxiaoyuan新虫 (初入文坛)
|
[求助]
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 ] |
回复此楼» 猜你喜欢
能否退出参与的面上项目解除限项
已经有18人回复
-
如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》
已经有3人回复
-
今日不放榜？网传国自然预计 8 月 27 日可查结果
已经有21人回复
-
估计是周四
已经有4人回复
-
我面上完蛋了
已经有9人回复
-
2026国自然函评费到账
已经有18人回复
-
建议基金发布提前给出明确的时间点
已经有15人回复
-
范进中举一文的中心思想
已经有7人回复
-
朋友圈看到的
已经有9人回复
-
让我中一个面上吧！
已经有16人回复
yuwuxiaoyuan
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 444.8
- 帖子: 50
- 在线: 13.4小时
- 虫号: 2382452
- 注册: 2013-03-27
- 专业: 工程热物理相关交叉领域
2楼2013-04-08 18:42:34
ftai
金虫 (著名写手)
- 应助: 19 (小学生)
- 金币: 7836.4
- 散金: 574
- 红花: 7
- 帖子: 1422
- 在线: 356.1小时
- 虫号: 2101685
- 注册: 2012-11-02
- 性别: GG
- 专业: 信息安全
3楼2013-04-08 19:06:49
yuwuxiaoyuan
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 444.8
- 帖子: 50
- 在线: 13.4小时
- 虫号: 2382452
- 注册: 2013-03-27
- 专业: 工程热物理相关交叉领域
4楼2013-04-08 21:14:51
ftai
金虫 (著名写手)
- 应助: 19 (小学生)
- 金币: 7836.4
- 散金: 574
- 红花: 7
- 帖子: 1422
- 在线: 356.1小时
- 虫号: 2101685
- 注册: 2012-11-02
- 性别: GG
- 专业: 信息安全
5楼2013-04-08 21:57:33
yuwuxiaoyuan
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 444.8
- 帖子: 50
- 在线: 13.4小时
- 虫号: 2382452
- 注册: 2013-03-27
- 专业: 工程热物理相关交叉领域
6楼2013-04-08 22:55:02
20