24小时热门版块排行榜

返回列表

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

birdcanfly

金虫 (正式写手)

应助: 8 (幼儿园)
金币: 735.7
散金: 406
红花: 3
帖子: 320
在线: 105.9小时
虫号: 399057
注册: 2007-06-11
性别: GG
专业: 生物医用高分子材料

[求助] 微分方程ode45求解，最小二乘法优化微分方程参数，程序运行求助已有1人参与

仿照下面的帖子（http://muchong.com/bbs/viewthread.php?tid=6425538&authorid=1122189）进行了程修改，求解微分方程组的系数，微分方程组中包含了指前系数和活化能。
现在运行程序，matlab跑起来没完，无结果，无报错。请求高手指点，协助调通程序，50金币相谢！
***********************
程序代码如下：

function k1k2k3
format long
clear all
clc

k0 = [1e40 1e20 1e40 1e30 100e3 50e3 100e3 100e3 0.8];
lb = [0 0 0 0 0 0 0 0 0];
ub = [inf inf inf inf inf inf inf inf 1];

data=...
[
0.00 100.000 0.00000 0.00000 0.00000
0.21 78.8686 18.802 1.99667 0
0.22 59.2346 37.9368 2.66223 0.332779
0.23 48.0865 46.7554 4.49251 0.166389
0.24 37.1048 57.0715 5.32446 0.665557
0.25 29.6173 61.2313 8.31947 0.831947
0.253 31.1148 59.5674 9.65058 0.831947
0.258 29.7837 61.2313 7.8203 0.998336
0.26 22.629 63.3943 12.3128 1.16473
0.263 24.792 60.3993 13.4775 1.33111
0.266 21.9634 62.396 13.9767 1.4975
0.27 21.2978 63.228 13.1448 1.83028
0.272 24.6256 59.0682 14.3095 1.66389
0.277 17.3045 58.4027 21.1314 3.1614
0.28 17.4709 58.7354 20.9651 2.82862
0.29 23.7937 49.9168 22.4626 3.32779
0.3 18.3028 43.7604 31.1148 6.48918
];
x0=data(1,2:end);
tspan=data(:,1)';
yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5)];

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.9f ± %.9f\n',k(4),ci(4,2)-k(4))
fprintf('\tEa1 = %.9f ± %.9f\n',k(5),ci(5,2)-k(5))
fprintf('\tEa2 = %.9f ± %.9f\n',k(6),ci(6,2)-k(6))
fprintf('\tEa3 = %.9f ± %.9f\n',k(7),ci(7,2)-k(7))
fprintf('\tEa4 = %.9f ± %.9f\n',k(8),ci(8,2)-k(8))
fprintf('\ta  = %.9f ± %.9f\n',k(9),ci(9,2)-k(9))
fprintf('  The sum of the squares is: %.9e\n\n',resnorm)

ts=0

(max(tspan)-min(tspan))/100):max(tspan);
[ts ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,3) data(:,4)];
I100=100*ones(16,1);
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo');
hold on
plot(ts,ys(:,2)+ys(:,3),'r',tspan,yy(:,2),'r*');
plot(ts,ys(:,4),'k',tspan,yy(:,3),'k+');
plot(ts,I100(:,1)-ys(:,1)-ys(:,2)-ys(:,3)-ys(:,4),'g',tspan,yy(:,4),'g<')
legend('C1的计算值','C1的实验值','C2的计算值*5','C2的实验值*5','C3的计算值*5','C3的实验值*5','C4的计算值','C4的实验值')

function f = ObjFunc(k,tspan,x0,yexp) % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);% 求解常微分方程，其中tspan为t的取值点，x0为微分方程组的初始值，k为微分方程的系数，返回t为微分方程组解的取值点，Xsim为微分方程组的解
Xsim1=Xsim(:,1);%提取Xsim的第一列
Xsim2=Xsim(:,2);%提取Xsim的第二列
Xsim3=Xsim(:,3);%提取Xsim的第三列
Xsim4=Xsim(:,4);%提取Xsim的第四列

ysim(:,1) = Xsim1(2:end);%微分方程的第一个变量，从第二个解到最后一个解赋值给ysim的第一列
ysim(:,2) = Xsim2(2:end);%微分方程的第二个变量，从第二个解到最后一个解赋值给ysim的第二列
ysim(:,3) = Xsim3(2:end);%微分方程的第三个变量，从第二个解到最后一个解赋值给ysim的第三列
ysim(:,4) = Xsim4(2:end);%微分方程的第四个变量，从第二个解到最后一个解赋值给ysim的第四列

I100=100*ones(16,1);
f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)+ysim(:,3)-yexp(:,2)) (ysim(:,4)-yexp(:,3)) (I100-ysim(:,1)-ysim(:,2)-ysim(:,3)-ysim(:,4)-yexp(:,4))];%形成目标优化函数

function dCdt = KineticsEqs(t,C,k) % ODE模型方程，C为浓度，k为微分方程中的系数，t为导数
a=k(9);
Ea1=k(5);
Ea2=k(6);
Ea3=k(7);
Ea4=k(8);
R=8.314;
Ttime=(-960.5*t^2+842.1*t+56.66+273.15);
dC1dt =-k(1)*exp(-Ea1/(R*Ttime))*C(1);%第一个微分方程，等号前面是C(1)对t的微分
dC2dt =k(1)*exp(-Ea1/(R*Ttime))*a*C(1)-k(2)*exp(-Ea2/(R*Ttime))*C(2);%第二个微分方程，等号前面是C(2)对t的微分
dC3dt =k(2)*exp(-Ea2/(R*Ttime))*C(2)-k(3)*exp(-Ea3/(R*Ttime))*C(3);%第三个微分方程，等号前面是C(3)对t的微分
dC4dt =k(3)*exp(-Ea3/(R*Ttime))*C(3)-k(4)*exp(-Ea4/(R*Ttime))*C(4);%第四个微分方程，等号前面是C(4)对t的微分
dCdt = [dC1dt; dC2dt;dC3dt;dC4dt];%微分方程组

回复此楼