版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

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];%微分方程组

回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

1楼 2014-09-18 14:11:06

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

whqs8426212

铜虫 (正式写手)

应助: 5 (幼儿园)
金币: 13.2
散金: 375
红花: 7
帖子: 782
在线: 685.7小时
虫号: 2572504
注册: 2013-07-30
专业: 数论

【答案】应助回帖

感谢参与，应助指数 +1

两个边界值范围缩小，试试

[ 发自小木虫客户端 ]

赞一下

回复此楼

5楼2014-09-20 18:18:39

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

查看全部 6 个回答

dikeway

新虫 (小有名气)

应助: 1 (幼儿园)
金币: 320.3
散金: 1
红花: 1
帖子: 102
在线: 34.4小时
虫号: 3414256
注册: 2014-09-13
性别: GG
专业: 经济地理学

帮顶一下，我的也没回复呢

赞一下

回复此楼

2楼2014-09-18 22:12:01

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

birdcanfly

金虫 (正式写手)

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

可以顺利运行，求出了参数，可以很好的拟合文献中的数据，但是参数和文献报道的参数差好几个数量级。
*********************
function k1k2k3
format short e
clear all
clc
k0 = [1.0 1.0 1.0 1.0 139 59 169 100 0.8];
lb = [0 0 0 0 0 0 0 0 0];
ub = [inf inf inf inf inf inf inf inf 1];

data=...
[
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 = %e ± %e\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %e ± %e\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %e ± %e\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %e ± %e\n',k(4),ci(4,2)-k(4))
fprintf('\tEa1 = %e ± %e\n',k(5),ci(5,2)-k(5))
fprintf('\tEa2 = %e ± %e\n',k(6),ci(6,2)-k(6))
fprintf('\tEa3 = %e ± %e\n',k(7),ci(7,2)-k(7))
fprintf('\tEa4 = %e ± %e\n',k(8),ci(8,2)-k(8))
fprintf('\ta = %f ± %f\n',k(9),ci(9,2)-k(9))

ts=0.210

(max(tspan)-min(tspan))/100):max(tspan);
[ts ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,3) data(:,4) data(:,5)];
I100=100*ones(101,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-ys(:,1)-ys(:,2)-ys(:,3)-ys(:,4),'g',tspan,yy(:,4),'g<')
legend('C1的计算值','C1的实验值','C2的计算值','C2的实验值','C3的计算值','C3的实验值','C4的计算值','C4的实验值')

function f = ObjFunc(k0,tspan,x0,yexp) % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k0);% 求解常微分方程，其中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(15,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];%微分方程组

赞一下

回复此楼

3楼2014-09-19 13:17:46

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

birdcanfly

金虫 (正式写手)

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

引用回帖:

3楼: Originally posted by birdcanfly at 2014-09-19 13:17:46
可以顺利运行，求出了参数，可以很好的拟合文献中的数据，但是参数和文献报道的参数差好几个数量级。
*********************
function k1k2k3
format short e
clear all
clc
k0 = ;
lb = ;
ub = ;

data= ...

将参数的初始值设定为文献中报道的数值，报错如下：

Local minimum possible.

lsqnonlin stopped because the size of the current step is less than
the default value of the step size tolerance.

<stopping criteria details>

Warning: Matrix is singular to working precision.
> In nlparci at 104
In ODEParmafit_non_test at 33

赞一下

回复此楼

4楼2014-09-19 13:22:20

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

查看全部 6 个回答

最具人气热帖推荐 [查看全部]		作者	回/看	最后发表

[考研] 材料学硕297已过四六级求调剂推荐 +6	adaie 2026-03-19	6/300	2026-03-20 14:12 by 星空星月
[考研] 295复试调剂 +7	简木ChuFront 2026-03-19	7/350	2026-03-20 14:11 by 星空星月
[考研] 本人考085602 化学工程专硕 +18	不知道叫什么！ 2026-03-15	20/1000	2026-03-20 13:52 by danney002
[考研] 材料专硕英一数二306 +6	z1z2z3879 2026-03-18	6/300	2026-03-20 08:49 by xingguangj
[考研] 321求调剂 +8	何润采123 2026-03-18	10/500	2026-03-19 16:46 by 何润采123
[考研] 一志愿北京化工大学0703化学318分，有科研经历，求调剂 +3	一瓶苯甲酸 2026-03-14	3/150	2026-03-19 15:17 by 尽舜尧1
[考研] 085600材料与化工调剂 324分 +10	llllkkkhh 2026-03-18	12/600	2026-03-19 14:33 by llllkkkhh
[考研] 324分 085600材料化工求调剂 +3	llllkkkhh 2026-03-18	3/150	2026-03-19 14:22 by houyaoxu
[考研] 一志愿福大288有机化学，求调剂 +3	小木虫200408204 2026-03-18	3/150	2026-03-19 13:31 by houyaoxu
[考研] 328求调剂，英语六级551，有科研经历 +3	生物工程调剂 2026-03-17	7/350	2026-03-18 20:41 by Wangjingyue
[考研] 材料专业求调剂 +5	hanamiko 2026-03-18	5/250	2026-03-18 20:19 by 星空星月
[考研] 070300化学319求调剂 +6	锦鲤0909 2026-03-17	6/300	2026-03-18 13:22 by Iveryant
[考研] 299求调剂 +5	△小透明* 2026-03-17	5/250	2026-03-18 11:49 by 尽舜尧1
[考研] 材料，纺织，生物（0856、0710），化学招生啦 +3	Eember. 2026-03-17	9/450	2026-03-18 10:28 by Eember.
[考研] 268求调剂 +6	简单点0 2026-03-17	6/300	2026-03-18 09:04 by 无际的草原
[考研] 301求调剂 +9	yy要上岸呀 2026-03-17	9/450	2026-03-18 08:58 by 无际的草原
[考研] 277调剂 +5	自由煎饼果子 2026-03-16	6/300	2026-03-17 19:26 by 李leezz
[考研] 283求调剂 +3	听风就是雨； 2026-03-16	3/150	2026-03-17 07:41 by 热情沙漠
[考研] 东南大学364求调剂 +5	JasonYuiui 2026-03-15	5/250	2026-03-16 21:28 by 木瓜膏
[考研] 0703 物理化学调剂 +3	我可以上岸的对� 2026-03-13	5/250	2026-03-16 10:50 by 我可以上岸的对�

24小时热门版块排行榜

birdcanfly

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

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

whqs8426212

【答案】应助回帖

dikeway

birdcanfly

birdcanfly

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