求助代謝動(dòng)力學(xué)系數(shù)模擬代碼 1stOpt或者M(jìn)ATLAB
求代謝動(dòng)力學(xué)系數(shù)模擬
一級(jí)代謝物的動(dòng)力學(xué)方程:dC2/dt=k1*C1-k2*C2
初始條件:t=0,C2=0
且C1=exp(-A*t)。A=0.2779
我嘗試了用1stOpt(破解版)和MATLAB ODE方法,都沒(méi)成功,想請(qǐng)教一下大神。
另外t不是嚴(yán)格的等差數(shù)列,取值如: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使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n')
fprintf('\t待擬合參數(shù) k1 = %.6f\n',k(1))
fprintf('\t待擬合參數(shù) 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決定系數(shù)R-Square = %.6f',R2);
figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('計(jì)算值','實(shí)驗(yàn)值','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的,沒(méi)有做修改)
%---------------------------------------------------------
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;
求求啦,被這個(gè)問(wèn)題卡了兩個(gè)多月了,不知道怎么解出k1 k2
返回小木蟲(chóng)查看更多
京公網(wǎng)安備 11010802022153號(hào)
參數(shù)擬合C2缺少數(shù)據(jù)
1stOpt容易實(shí)現(xiàn),1.5不支持微分方程擬合,需要下載5.0版本的
謝謝回復(fù)。
有C2的數(shù)據(jù),因?yàn)镃1服從指數(shù)方程,直接把方程寫(xiě)在程序里了。C2的數(shù)據(jù)就寫(xiě)在了下面的數(shù)據(jù)表里。
能分享5.0版的1stOpt下載么?我在網(wǎng)上實(shí)在找不到了,我這里也沒(méi)有校內(nèi)bbs之類(lèi)的東西。。
謝謝了
C2數(shù)據(jù)給出來(lái)看看
不是數(shù)據(jù)的問(wèn)題,是程序不能運(yùn)行的問(wèn)題。點(diǎn)運(yùn)行之后就轉(zhuǎn)圈,輸出那里也沒(méi)反應(yīng)。
將C 的數(shù)據(jù)補(bǔ)上就可以運(yùn)行了。
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
好像不對(duì)頭。
大哥加個(gè)qq啊,我好請(qǐng)教你?br>,