24小时热门版块排行榜    

查看: 2651  |  回复: 10
【悬赏金币】回答本帖问题,作者Evoly_z将赠送您 88 个金币
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

Evoly_z

银虫 (小有名气)

[求助] 求助代谢动力学系数模拟代码 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
回复此楼
RELEASE/ME/FROM/THE/PRISON!
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

独孤神宇

版主 (知名作家)

【答案】应助回帖

将C 的数据补上就可以运行了。

Parameters k1[0,100], k2[0,100];
Variable t, C;
InitialODEValue t=0,c=0;
ODEFunction c'=k1*exp(-0.2779*t)-k2*c;
Data;
//t    C
   0    xxx
   1    xxx
   2    xxx
   4    xxx
   6    xxx
   10  xxx
   15  xxx
   24  xxx
数值计算
6楼2021-03-25 09:23:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 11 个回答

独孤神宇

版主 (知名作家)

【答案】应助回帖

感谢参与,应助指数 +1
参数拟合C2缺少数据
1stOpt容易实现,1.5不支持微分方程拟合,需要下载5.0版本的

发自小木虫Android客户端
数值计算
2楼2021-03-24 07:20:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Evoly_z

银虫 (小有名气)

引用回帖:
2楼: Originally posted by 独孤神宇 at 2021-03-24 07:20:06
参数拟合C2缺少数据
1stOpt容易实现,1.5不支持微分方程拟合,需要下载5.0版本的

谢谢回复。

有C2的数据,因为C1服从指数方程,直接把方程写在程序里了。C2的数据就写在了下面的数据表里。

能分享5.0版的1stOpt下载么?我在网上实在找不到了,我这里也没有校内bbs之类的东西。。

谢谢了
RELEASE/ME/FROM/THE/PRISON!
3楼2021-03-24 11:13:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
C2数据给出来看看
4楼2021-03-24 16:12:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
不应助 确定回帖应助 (注意:应助才可能被奖励,但不允许灌水,必须填写15个字符以上)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 应该是93bebmhtak前后十一个字符比较关键 +22 Lanmanbaby 2026-08-09 36/1800 2026-08-12 22:50 by sdfapple719
[基金申请] 关于Filecode分析方法 +8 majunge000 2026-08-10 9/450 2026-08-12 22:17 by 虫友是什么虫
[基金申请] 2019年青年基金涵评意见,大家看看几个A,几个B? +10 Tide man 2026-08-11 10/500 2026-08-12 20:59 by 陈晨晨陈啊
[基金申请] 好奇怪的filecode +5 布布和一二 2026-08-08 6/300 2026-08-12 16:11 by 云上清扬
[基金申请] 长年满屏的广告,版主太不责任了。基金也等的急 +4 gltch 2026-08-12 4/200 2026-08-12 16:08 by eulota
[论文投稿] 职称评审,求友友推荐见刊最快的期刊 +4 工厂打螺丝 2026-08-08 4/200 2026-08-12 09:18 by hansi2025
[基金申请] 综述论文作为代表作会不会影响评审专家的印象分? +11 yufeiwaner 2026-08-09 13/650 2026-08-12 08:17 by yufeiwaner
[基金申请] 小木虫上这么多卖论文的,真有人买论文么?感觉没必要啊 +8 Tide man 2026-08-10 9/450 2026-08-11 20:39 by beefly
[基金申请] 有时候,自然基金真的不能太认真 (我的申报经验) +3 majunge000 2026-08-11 4/200 2026-08-11 20:13 by lch2012
[基金申请] 帮忙看看fileCode +7 wwncly 2026-08-10 13/650 2026-08-11 19:36 by 冰心玉壶晴
[基金申请] 为什么网上很多人说本周 12号出结果 +6 瞬息宇宙 2026-08-10 7/350 2026-08-11 19:25 by Tide man
[基金申请] 什么时候出结果,有咨询渠道??? +3 Tide man 2026-08-11 3/150 2026-08-11 17:54 by kudofaye
[基金申请] 关于代码变化问题,想知道的进来 +17 且听虎啸 2026-08-07 24/1200 2026-08-10 18:35 by zhangduo2008
[基金申请] 国自然结果 +4 Vierhys 2026-08-10 8/400 2026-08-10 15:06 by Vierhys
[基金申请] 2026国自然放榜时间 +9 布布和一二 2026-08-08 9/450 2026-08-10 11:22 by xxxx2020
[基金申请] 这样的filecode谁见过 +11 布布和一二 2026-08-08 22/1100 2026-08-10 11:10 by wmfsnow
[基金申请] filecode与中标关系的预测 +5 布布和一二 2026-08-07 5/250 2026-08-09 16:15 by 袁向阳007
[基金申请] 关于filecode,很负责任的告诉大家 +6 爱看书的可乐 2026-08-08 7/350 2026-08-08 22:13 by a_niu
[基金申请] 关于filecode +4 布布和一二 2026-08-07 7/350 2026-08-07 22:55 by zhanghaozhu
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
信息提示
请填处理意见