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

登录注册

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

24小时热门版块排行榜

北京石油化工学院2026年研究生招生接收调剂公告

返回列表

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

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

回复此楼

» 猜你喜欢

【新加坡】纳米电子器件项目组有“联合培养博士生”名额已经有2人回复
重庆交大材料学院26年硕士招生复试实施细则已出。欢迎加入光子学微结构与器件课题组。已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有279人回复
北工大宋晓艳教授团队招收2026年9月入学博士【学博和专博名额均有】！已经有1人回复
天津城建大学理学院化学、环境工程、材料物理等方向接收10人已经有0人回复
天津城建大学理学院化学、环境工程、材料物理等方向接收调剂，10人已经有0人回复
计算与实验科学工程国际会议(8月7-11日,香港) 二维材料分会S5与声学换能器分会S8 已经有0人回复
0702一志愿吉大B区求调剂有论文已经有4人回复
欢迎调剂到赣南师范大学智能制造与未来能源学院已经有10人回复
0702一志愿吉大B区求调剂本科期间发表一篇Sci 已经有4人回复

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

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

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

birdcanfly

金虫 (正式写手)

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

引用回帖:

5楼: Originally posted by whqs8426212 at 2014-09-20 18:18:39
两个边界值范围缩小，试试

如果不知道参数值，如何缩小边界？刚才试了试退火算法，也是一样模拟值和给出的初始值关系密切。

赞一下

回复此楼

6楼2014-09-23 18:08:55

已阅回复此楼关注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 个回答

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

[考研] 一志愿中海洋320化学工程与技术学硕求调剂 +3	披星河 2026-03-30	3/150	2026-03-30 23:38 by 果果妈咪
[考研] 求调剂，一志愿南京航空航天大学，080500材料科学与工程学硕，总分289分 +9	@taotao 2026-03-29	9/450	2026-03-30 22:29 by 我是小康
[考研] 085701环境工程求调剂 +11	多久上课 2026-03-27	12/600	2026-03-30 21:21 by 研究僧导导
[考研] 一志愿北京化工大学材料与化工（085600）296求调剂 +25	稻妻小编 2026-03-26	25/1250	2026-03-30 20:11 by 滴滴上岸呀
[考研] 367求调剂 +5	芋泥啵啵… 2026-03-28	5/250	2026-03-30 19:56 by 无际的草原
[考研] 材料学硕333求调剂 +14	北道巷 2026-03-24	14/700	2026-03-30 18:59 by 源_2020
[考研] 求调剂 +10	家佳佳佳佳佳 2026-03-29	10/500	2026-03-30 18:34 by 544594351
[考研] 求调剂 +7	研研，接电话 2026-03-24	8/400	2026-03-30 14:14 by 白云朵朵飞
[考研] 085602 化学工程专硕 340分求调剂 +4	qianbai11 2026-03-29	4/200	2026-03-30 11:34 by 唐沐儿
[考研] 求化学调剂 +11	wulanna 2026-03-28	11/550	2026-03-30 10:59 by 探123
[考研] 337求调剂 +6	《树》 2026-03-29	6/300	2026-03-30 10:15 by herarysara
[考研] 275求调剂 +15	Micky11223 2026-03-25	20/1000	2026-03-29 20:44 by 唐沐儿
[考研] 085600，专业课化工原理，321分求调剂 +5	大馋小子 2026-03-28	5/250	2026-03-29 08:56 by qingfeng258
[考研] 081200-11408-276学硕求调剂 +6	崔wj 2026-03-26	6/300	2026-03-29 01:11 by hanserlol
[考研] 药学105500求调剂 +3	Ssun。。 2026-03-28	3/150	2026-03-28 11:24 by lxf170613
[考研] 一志愿南师大0703化学 275求调剂 +4	Ripcord上岸 2026-03-27	4/200	2026-03-27 17:00 by zhyzzh
[考研] 305求调剂 +5	哇卢卡库 2026-03-26	5/250	2026-03-27 14:01 by laoshidan
[考研] 085602化学工程求调剂。 +4	平乐乐乐 2026-03-26	4/200	2026-03-26 17:57 by fmesaito
[考研] 340求调剂 +5	话梅糖111 2026-03-24	5/250	2026-03-25 06:53 by ilovexiaobin
[考研] 080500求调剂 +3	zzzzfan 2026-03-24	3/150	2026-03-24 16:38 by barlinike

24小时热门版块排行榜

birdcanfly

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

» 猜你喜欢

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

birdcanfly

dikeway

birdcanfly

birdcanfly

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