![]() |
[求助]这是从一篇文章摘下的程序,我是一名maltlab初学者,应急需使用此程序,恳请帮助!
这有两个程序,上面的可以执行,下面的却无法运行,请大牛们帮忙看看,先谢了!
第一个程序 clear all clc global Cam t=[0 4 8 12 16 20 24 28 32 36 40 44 48 52]; CAm=[0.006 0.047 0.511 2.070 3.170 3.944 4.283 4.236 4.163 3.965 3.804 3.681 3.645 3.585]'; %动力学数据 % 非线性拟合 beta0 = [0.5816 4.250]; tspan = [0 4 8 12 16 20 24 28 32 36 40 44 48 52]; CA0 = 0.006; lb=[0.2 0.5]; ub=[0.6 4.0] [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@OptObjFunc5,beta0,lb,ub,[],tspan,CA0,CAm) ci = nlparci(beta,resid,jacobian) % 拟合效果图(实验与拟合的比较) [t4plot CA4plot] = ode45(@KineticsEqs5,[tspan(1) tspan(end)],CA0,[],beta); plot(t,CAm,'bo',t4plot,CA4plot,'k-') legend('Exp','Model') xlabel('时间t, s') ylabel('浓度C_A, mol/L') % 残差关于拟合值的残差图 [t CAc] = ode45(@KineticsEqs5,tspan,CA0,[],beta); figure plot(CAc,resid,'*') xlabel('浓度拟合值(mol/L)') ylabel('残差R (mol/L)') refline(0,0) % 参数辨识结果 fprintf('Estimated Parameters:\n') fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1)) fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2)) function f = OptObjFunc5(beta,tspan,CA0,CAm) [t CAc] = ode45(@KineticsEqs5,tspan,CA0,[],beta); f = CAc - CAm; % ------------------------------------------------------------------ function dCAdt = KineticsEqs5(t,CA,beta) dCAdt = beta(1)*(1-CA/beta(2))*CA; % k= beta(1), n= beta(2) 第二个程序 k0=[0.5818 20.0957 0.028907]; %μmxa、Yx/s、beta参数初值 x0=[0.006 36.001 0.006]; %x0,s0,p0 t=[0 4 8 12 16 20 24 28 32 36 40 44 48 52]; tspan=[0 4 8 12 16 20 24 28 32 36 40 44 48 52]; %yexp实验数据[xl x2 x3] yexp=[[0.006 0.047 0.511 2.070 3.170 3.944 4.283 4.236 4.1633 .965 3.804 3.681 3.645 3.585]; [36.001 35.788 34.735 32.215 27.187 20.593 10.802 5.850 3.129 2.251 1.850 1.499 1.111 0.960]; [0.006 0.007 0.008 0.009 0.030 0.139 0.605 1.033 1.621 2.251 2.544 2.840 3.249 3.300]]’: lb=[0.2 0.05 0.01]; ub=[0.6 1 0.031]; %参数下、上限 % [k, resnorm,resid,exitflag,output,lambda, jacobian]=… lsnqonlin(@ObjFunc,k0,lb,ub,[],x0,yexp) ci=nlparci(k, resid, jacoban) % yl=[0.006 0.047 0.511 2.070 3.170 3.944 4.283 4.236 4.1633 .965 3.804 3.681 3.645 3.585]'; y2=[36.001 35.788 34.735 32.215 27.187 20.593 10.802 5.850 3.129 2.251 1.850 1.499 1.111 0.960]'; y3=[0.006 0.007 0.008 0.009 0.030 0.139 0.605 1.033 1.621 2.251 2.544 2.840 3.249 3.300]'; [t4plot x4plot]=ode45(@kineticseqs,[tspan(1) tspan(end)],x0,[],k); plot(t1,yl,'bo', tl,y2', 'g*',tl,y3,'r*', t4plot,x4Plot,'k-'), legend('Exp-biomass','Exp-Na-citrate','Exp-产物活性单位','Model'), xlabel('time(h)'),ylabel('菌体浓度,Na-citrate(g/l),产物活性单位(u)') % fprintf('Estimated Parameter\n') fprintf('\tk=%.4f±%.4f\n', k(1), ci(1,2)-k(1)) fprintf('\tk=%.4f±%.4f\n', k(2), ci(2,2)-k(2)) fprintf('\tk=%.4f±%.4f\n', k(1), ci(1,2)-k(3)) % [t x]=ode45(@kineticseqs,tspan,x0,[],k); figure,plot(x,resid,'*') xlabel('菌体、Na-citrate浓度拟合值(g/l)、产物活性单位(U)'),ylbael('残差R(g/l)') M文件: function f=ObjFunc (k,x0,yexp) tspan=[0 4 8 12 16 20 24 28 32 36 40 44 48 52]'; [t x]=ode45(@kineticseqs,tspan,x0,[],k); y(:,1)=x(:,l); y(:,2)=x(:,2); y(:,3)=x(:,3); %? fl=y(:,l)-yexp(:,1); f2=y(:,2)-yexp(:,2); f3=y(:,3)-yexp(:,3); f=[f1, f2, f3]; M文件: function dxdt=kineticseqs (t,x,k) %模型方程 dxdt=[k(l)*x(l)*(1.0-x(l)/4.25044) -k(l)/k(2)*x(l)*(1.0-x(l)/4.25044) k(3)*x(l)]; |
回复: [求助]这是从一篇文章摘下的程序,我是一名maltlab初学者,应急需使用此程序,恳请帮助!
你这个文章是什么?能传一下么?
|
所有时间均为北京时间。现在的时间是 13:13。 |
Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.