24小时热门版块排行榜

返回列表

【悬赏金币】回答本帖问题，作者Evoly_z将赠送您 88 个金币

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

Evoly_z

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 109.5
帖子: 51
在线: 30.9小时
虫号: 2116918
注册: 2012-11-09
性别: GG
专业: 环境污染化学

[求助] 求助代谢动力学系数模拟代码 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人回复
求取一些关于纳米材料和纳米技术相关的英文PPT。已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有298人回复
【复旦大学】二维材料方向招收2026年博士研究生1名已经有0人回复
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组2026级硕/博/博后招生已经有10人回复
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生已经有19人回复
反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
26申博推荐：南京航空航天大学国际前沿科学研究院光学方向招收博士生！已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复

RELEASE/ME/FROM/THE/PRISON!

1楼2021-03-24 02:33:08

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wlfc

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 77
帖子: 19
在线: 12小时
虫号: 26920802
注册: 2021-08-01

【答案】应助回帖

用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