| 查看: 2674 | 回复: 10 | |||
| 【悬赏金币】回答本帖问题，作者Evoly_z将赠送您 88 个金币 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[求助]
求助代谢动力学系数模拟代码 1stOpt或者MATLAB 已有3人参与
|
|||
 求代谢动力学系数模拟一级代谢物的动力学方程：dC2/dt=k1*C1-k2*C2 初始条件：t=0，C2=0 且C1=exp（-A*t）。A=0.2779 我尝试了用1stOpt(破解版)和MATLAB ODE方法，都没成功，想请教一下大神。 另外t不是严格的等差数列，取值如：t=0，1，2，4，6，10，15，24 1stOpt代码： Title Kinetic_ave Parameters k1[0,100], k2[0,100]; Variable t, C; StartProgram var i:integer; begin for i:=0 to DataLength -1 do begin if i ==0 C=0; else C:=C[i-1]+k1*(t-t[i-1])*exp(-0.2779 *t) - k2*C*(t-t[i-1]);end; EndProgram; Data; //t C 0 xxx 1 xxx 2 xxx 4 xxx 6 xxx 10 xxx 15 xxx 24 xxx Matlab代码： function ODE_ave clear all;clc format long aveall; t=T_h(  ;yexp=OLEave( ;k0=[1 1]; y0=0; lb=[0 0]; ub=[+inf +inf]; yy=[y0 yexp']; tspan=0:1:24; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\t待拟合参数 k1 = %.6f\n',k(1)) fprintf('\t待拟合参数 k2 = %.6f\n',k(2)) fprintf(' \t残差平方和= %.6f\n\n',resnorm) ts=0:1:24; [ts ys]=ode45(@KineticsEqs,ts,y0,[],k); [ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k); y=XXsim(2:end); xexp=yexp; R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2); fprintf('\n\t决定系数R-Square = %.6f',R2); figure(1) plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best'); yr=y-yexp; figure(2) plot(tspan(2:end),yr,'r*',[-1 15],[0 0]),axis([-1 15 -0.5 0.5]); figure(3) plot(yexp,y,'ro',[21 29],[21 29],'b-'); (作图这块儿是copy的，没有做修改) %--------------------------------------------------------- function f = ObjFunc(k,tspan,y0,yexp) [t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ; ysim = Xsim(2:end); size(ysim); size(yexp); f=ysim(1,1)+ysim(2,1)+ysim(4,1)+ysim(6,1)+ysim(10,1)+ysim(15,1)+ysim(24,1) - sum(yexp(:,1)); %---------------------------------------------------------- function dydt = KineticsEqs(t,y,k) beta(1)=k(1); beta(2)=k(2); dydt = beta(1)*exp(-0.2779*t)-beta(2)*y; 求求啦，被这个问题卡了两个多月了，不知道怎么解出k1 k2 |
回复此楼» 猜你喜欢
如何从铁电相到顺电相。
已经有1人回复
-
各位大佬，求一份最新Wien2k版本作为学习用途
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有224人回复
-
南科大活性流体和软物质课题组诚招2027级博士、硕士研究生和博后
已经有42人回复
-
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有130人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
PVA溶解
已经有2人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
【答案】应助回帖
|
用OpenLu求解，Lu脚本代码： !!!using["luopt","math"]; //使用命名空间 f(x,y,dy, params::k1,k2)= { dy=k1*(0.701*exp(-3.211*x)+0.299*exp(-0.067*x))-k2*y, 0 //必须返回0 }; 目标函数(_k1,_k2 : i,s,tyz : tyArray,tA,max,k1,k2)= { k1=_k1, k2=_k2, //传递优化变量 //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度 tyz=gsl_ode[@f,nil,0.0,tA,ra1(0), 1e-6, 1e-6, 2, 1e-6,50], i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2}, s }; main(::tyArray,tA,max)= { tyArray=matrix{ //存放实验数据xi,yi "0 0 0.33 0.061043 1 0.03675 2 0.05932 4 0.095993 6 0.072057 10 0.05085 15 0.04678 24 0.047673 30 0.034973 48 0.030375 72 0" }, len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数，tA取矩阵的列 Opt1[@目标函数] //Opt1函数全局优化 }; 结果： 0.1325616120497096 0.358508535326814 3.242596823546582e-003 |
11楼2021-08-09 21:27:36
独孤神宇
版主 (知名作家)
- 应助: 490 (硕士)
- 贵宾: 0.008
- 金币: 31016.3
- 散金: 802
- 红花: 122
- 沙发: 1
- 帖子: 5600
- 在线: 856.5小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
2楼2021-03-24 07:20:06
3楼2021-03-24 11:13:13
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
4楼2021-03-24 16:12:00
25