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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 微分方程ode45求解，最小二乘法优化微分方程参数，程序运行求助
51/1返回列表
查看: 1858  |  回复: 5
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

birdcanfly

金虫 (正式写手)

[求助] 微分方程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的回帖

birdcanfly

金虫 (正式写手)
引用回帖:
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 个回答

dikeway

新虫 (小有名气)
帮顶一下，我的也没回复呢
赞一下 回复此楼
2楼2014-09-18 22:12:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

birdcanfly

金虫 (正式写手)
可以顺利运行，求出了参数，可以很好的拟合文献中的数据，但是参数和文献报道的参数差好几个数量级。
*********************
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的回帖

whqs8426212

铜虫 (正式写手)

【答案】应助回帖

感谢参与，应助指数 +1
两个边界值范围缩小，试试

[ 发自小木虫客户端 ]
赞一下 回复此楼
5楼2014-09-20 18:18:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 332求调剂 +4 ydfyh 2026-03-17 4/200 2026-03-21 02:20 by JourneyLucky
[考研] 085700资源与环境308求调剂 +12 墨墨漠 2026-03-18 13/650 2026-03-21 01:42 by JourneyLucky
[考研] 324分 085600材料化工求调剂 +4 llllkkkhh 2026-03-18 4/200 2026-03-21 01:24 by JourneyLucky
[考研] 材料专业求调剂 +6 hanamiko 2026-03-18 6/300 2026-03-21 00:24 by JourneyLucky
[考研] 296求调剂 +6 www_q 2026-03-18 10/500 2026-03-20 23:56 by JourneyLucky
[考研] 材料学硕297已过四六级求调剂推荐 +11 adaie 2026-03-19 11/550 2026-03-20 21:30 by laoshidan
[考研] 260求调剂 +3 朱芷琳 2026-03-20 3/150 2026-03-20 20:35 by 学员8dgXkO
[考研] 353求调剂 +3 拉钩不许变 2026-03-20 3/150 2026-03-20 19:56 by JourneyLucky
[考研] 环境工程调剂 +9 大可digkids 2026-03-16 9/450 2026-03-20 17:38 by 醉在风里
[考研] 085410人工智能专硕317求调剂（0854都可以） +4 xbxudjdn 2026-03-18 4/200 2026-03-20 09:07 by 不168
[考研] 求调剂 +3 暗涌afhb 2026-03-16 3/150 2026-03-20 00:28 by 河南大学校友
[考研] 288求调剂，一志愿华南理工大学071005 +5 ioodiiij 2026-03-17 5/250 2026-03-19 18:22 by zcl123
[考研] 266求调剂 +5 阳阳哇塞 2026-03-14 10/500 2026-03-19 15:08 by 阳阳哇塞
[考研] 0817调剂 +3 没有答案_ 2026-03-14 3/150 2026-03-19 09:51 by Xu de nuo
[考研] 311求调剂 +6 26研0 2026-03-15 6/300 2026-03-18 14:43 by haxia
[考研] 293求调剂 +11 zjl的号 2026-03-16 16/800 2026-03-18 08:10 by zhukairuo
[考博] 26博士申请 +3 1042136743 2026-03-17 3/150 2026-03-17 23:30 by 轻松不少随
[考研] 302求调剂 +4 小贾同学123 2026-03-15 8/400 2026-03-17 10:33 by 小贾同学123
[考研] 283求调剂 +3 听风就是雨； 2026-03-16 3/150 2026-03-17 07:41 by 热情沙漠
[考研] [导师推荐]西南科技大学国防/材料导师推荐 +3 尖角小荷 2026-03-16 6/300 2026-03-16 23:21 by 尖角小荷
信息提示
关闭
请填处理意见
关闭