| 查看: 1301 | 回复: 6 | ||
[求助]
求1stopt指点 已有1人参与
|
|
下面是我利用低版本的1stopt编写四阶龙格库塔法计算动力学方程，但是程序或许有问题我没有算出结果，请指点； Title "Kinetics"; //Parameters k1,k2,k3,k4,k5,k6,k7,k8,k9; //Variable t,F0,F,yH2,yDMM,yMe,yH2in,yDMMin,yMein; //Variable y1[output],y2[output],y3[output]; StartProgram [VB]; Dim iter As Integer dim yHPMin As Double dim yPDOin As Double dim yNPAin As Double dim F0 As Double dim F As Double Dim RK_i As Integer Dim h As Double Dim i As Integer dim j As Integer dim c As Double dim r1 As Double dim r2 As Double dim r3 As Double dim n_H2 As Double dim n_DMM As Double dim n_Me As Double dim y_H2 As Double dim y_DMM As Double dim y_Me As Double dim y_HPM As Double dim y_PDO As Double dim y_NPA As Double dim k1 As Double dim k2 As Double dim k3 As Double dim k4 As Double dim k5 As Double dim k6 As Double dim k7 As Double dim k8 As Double dim k9 As Double dim Y_Y01 As Double dim Y_Y02 As Double dim Y_Y03 As Double dim Y_Y04 As Double dim Y_Y05 As Double dim Y_Y06 As Double dim y_y1 As Double dim y_y2 As Double dim y_y3 As Double dim x_x1 As Double dim x_x2 As Double dim x_x3 As Double dim d_d1 As Double dim d_d2 As Double dim d_d3 As Double dim e_e1 As Double dim e_e2 As Double dim e_e3 As Double dim e_e4 As Double dim e_e5 As Double for iter to Datalenth-1 yHPMin=0 yPDOin=0 yNPAin=0 h=1.0008/20 Y_Y01=yHPMin Y_Y02=yPDOin Y_Y03=yNPAin Y_Y04=yH2in(iter) Y_Y05=yDMMin(inter) Y_Y06=yMein(inter) For RK_i=1 to 20 e_e1=0.5 e_e2=0.5 e_e3=1 e_e4=1 e_e5=0.5 x_x1=Y_Y01 x_x2=Y_Y02 x_x3=Y_Y03 For j=1 To 4 y_HPM=x_x1 y_PDO=x_x2 y_NPA=x_x3 n_H2=F0*yH2- 2*F*(y_HPM+y_PDO+y_NPA ) n_DMM=F0*yDMM- F*(y_HPM+y_PDO+y_NPA ) n_Me=F0*yMe+F*(y_HPM+y_PDO) y_H2= n_H2/F y_DMM= n_DMM/F y_ME= n_Me/F r1 = k7*k1*k2*y_H2*y_DMM /(1+k1*y_H2+k2*y_DMM+k3*y_Me+k4*y_HPM+k5*y_PDO+k6*y_NPA)^2 r2 = k8*k1*k4*y_H2*y_HPM/((1+k1*y_H2+k2*y_DMM+k3*y_Me+k4*y_HPM+k5*y_PDO+k6*y_NPA)^2 r3 = k9*k1*k5*y_H2*y_PDO/((1+k1*y_H2+k2*y_DMM+k3*y_Me+k4*y_HPM+k5*y_PDO+k6*y_NPA)^2 d_d1 = r1-r2 d_d2= r2-r3 d_d3 = r3 c=h*d_d1 if j=1 then x_x1=Y_Y01+e_e1*c y_y1=y_y1+e_e2*c/3 elseif j=2 then x_x1=Y_Y01+e_e2*c y_y1=y_y1+e_e3*c/3 elseif j=3 then x_x1=Y_Y01+e_e3*c y_y1=y_y1+e_e4*c/3 elseif j=4 then x_x1=Y_Y01+e_e4*c y_y1=y_y1+e_e5*c/3 end if c=h*d_d2 if j=1 then x_x2=Y_Y02+e_e1*c y_y2=y_y2+e_e2*c/3 elseif j=2 then x_x2=Y_Y02+e_e2*c y_y2=y_y2+e_e3*c/3 elseif j=3 then x_x2=Y_Y02+e_e3*c y_y2=y_y2+e_e4*c/3 elseif j=4 then x_x2=Y_Y02+e_e4*c y_y2=y_y2+e_e5*c/3 end if c=h*d_d3 if j=1 then x_x3=Y_Y03+e_e1*c y_y3=y_y3+e_e2*c/3 elseif j=2 then x_x3=Y_Y03+e_e2*c y_y3=y_y3+e_e3*c/3 elseif j=3 then x_x3=Y_Y03+e_e3*c y_y3=y_y3+e_e4*c/3 elseif j=4 then x_x3=Y_Y03+e_e4*c y_y3=y_y3+e_e5*c/3 end if Next Y_Y01=y_y1 Y_Y02=y_y2 Y_Y03=y_y3 Next y1(inter)=y_y1 y2(inter)=y_y2 y3(inter)=y_y3 Next EndProgram; Data; 2.05 0.48911 0.48786 0.92932 0.00258 0.06810 0.000031 0.00116 0.001137 1.37 0.73291 0.73179 0.93029 0.00155 0.06817 0.000078 0.00135 0.000573 1.62 0.61922 0.61815 0.91758 0.00173 0.08068 0.000104 0.00173 0.000475 1.98 0.50545 0.50453 0.89928 0.00187 0.09884 0.000289 0.00236 0.000444 1.21 0.82549 0.82425 0.91772 0.00158 0.08070 0.000414 0.00190 0.000280 1.49 0.67418 0.67281 0.89895 0.00224 0.09881 0.000641 0.00202 0.000292 2.70 0.37162 0.36991 0.81544 0.00531 0.17925 0.001304 0.00232 0.000368 1.02 0.97819 0.97574 0.92936 0.00254 0.06810 0.000206 0.00178 0.000361 0.82 1.22297 1.21986 0.92918 0.00273 0.06809 0.000465 0.00161 0.000189 0.97 1.03243 1.03056 0.91722 0.00213 0.08065 0.000790 0.00138 0.000149 1.19 0.84404 0.84146 0.89756 0.00379 0.09865 0.000972 0.00173 0.000210 2.16 0.46514 0.46287 0.81436 0.00663 0.17902 0.001744 0.00212 0.000349 |
回复此楼» 猜你喜欢
法国博士后职位
已经有0人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有178人回复
-
**
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有8人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有9人回复
-
桂林理工大学物理学专业招收调剂
已经有11人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- Table Curve 2D 数据拟合求助！
已经有5人回复
- 看了下1stOpt的介绍，感觉略显浮夸
已经有14人回复
- 求大哥大姐帮跑一下1stopt~新手啊，版本低啊
已经有6人回复
- Matlab用fsolve求解非线性方程，卡了三四天了跪求帮助！
已经有35人回复
- 高手指点：如何使用matlab使数据行与行间空一行
已经有6人回复
- 又会用proe5.0的吗 求大神 求教程
已经有2人回复
- 1stopt解方程组
已经有12人回复
- 在用origin拟合的过程中，程序编写这一块实在不会，求指教！
已经有10人回复
- 求高手帮忙拟合出方程
已经有7人回复
- 1stopt能进行隐函数非线性拟合吗？谁能帮拟合一下新型液体的粘度混合计算公式
已经有5人回复
- 求助如果用SPSS对已知模型进行非线形回归？
已经有7人回复
- 求大哥大姐帮跑一下1stopt~新手啊，版本低啊
已经有9人回复
- 用1stopt或者origin非线性拟合磁导率
已经有5人回复
- SPSS求非线性方程参数
已经有7人回复
- 1stopt4阶非线性常微分方程，帮小弟看一下，跪求。
已经有7人回复
- 1stopt或matlab如何求解以下的非线性方程并拟合出相应曲线？
已经有13人回复
- 麻烦高人指点一下1stopt和spss的非线性回归中weibull的两参数的公式应该怎么表达啊？
已经有9人回复
- 求SWMM软件高手指导
已经有15人回复
NicoleLas
木虫 (正式写手)
- 应助: 66 (初中生)
- 金币: 4288.7
- 红花: 26
- 帖子: 514
- 在线: 177.2小时
- 虫号: 2872321
- 注册: 2013-12-15
- 性别: GG
- 专业: 计算机软件
2楼2015-01-06 18:57:32
3楼2015-01-06 19:11:33
NicoleLas
木虫 (正式写手)
- 应助: 66 (初中生)
- 金币: 4288.7
- 红花: 26
- 帖子: 514
- 在线: 177.2小时
- 虫号: 2872321
- 注册: 2013-12-15
- 性别: GG
- 专业: 计算机软件
4楼2015-01-06 19:19:39
5楼2015-01-06 19:37:53
NicoleLas
木虫 (正式写手)
- 应助: 66 (初中生)
- 金币: 4288.7
- 红花: 26
- 帖子: 514
- 在线: 177.2小时
- 虫号: 2872321
- 注册: 2013-12-15
- 性别: GG
- 专业: 计算机软件
6楼2015-01-06 19:47:53
7楼2015-01-06 20:59:59
