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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 微分方程ode45求解，最小二乘法优化微分方程参数，程序运行求助
51/1返回列表
查看: 1852  |  回复: 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 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 317求调剂 +5 申子申申 2026-03-19 10/500 2026-03-20 15:58 by 蔡诚
[考研] 281求调剂（0805） +14 烟汐忆海 2026-03-16 25/1250 2026-03-20 15:47 by yuncha
[考研] 298-一志愿中国农业大学-求调剂 +9 手机用户 2026-03-17 9/450 2026-03-20 14:24 by 无懈可击111
[考研] 求调剂 +4 Mqqqqqq 2026-03-19 4/200 2026-03-20 14:15 by 星空星月
[考博] 招收博士1-2人 +3 QGZDSYS 2026-03-18 3/150 2026-03-20 11:58 by 呱呱呱呱叫
[考研] 材料080500调剂求收留 +6 一颗meteor 2026-03-13 6/300 2026-03-20 10:41 by EBSD
[考研] 求调剂 +3 暗涌afhb 2026-03-16 3/150 2026-03-20 00:28 by 河南大学校友
[论文投稿] 申请回稿延期一个月，编辑同意了。但系统上的时间没变，给编辑又写邮件了，没回复 10+3 wangf9518 2026-03-17 4/200 2026-03-19 23:55 by babero
[考研] 085600材料与化工调剂 324分 +10 llllkkkhh 2026-03-18 12/600 2026-03-19 14:33 by llllkkkhh
[考研] 材料工程专硕调剂 +5 204818@lcx 2026-03-17 6/300 2026-03-18 22:55 by 204818@lcx
[考研] 化学工程321分求调剂 +15 大米饭！ 2026-03-15 18/900 2026-03-18 14:52 by haxia
[考研] 312求调剂 +8 陌宸希 2026-03-16 9/450 2026-03-18 12:39 by Linda Hu
[考研] 0854，计算机类招收调剂 +3 胡辣汤放糖 2026-03-15 6/300 2026-03-18 12:09 by 上岸上岸……..
[考研] 0703化学336分求调剂 +6 zbzihdhd 2026-03-15 7/350 2026-03-18 09:53 by zhukairuo
[考研] 考研求调剂 +3 橘颂. 2026-03-17 4/200 2026-03-17 21:43 by 有只狸奴
[考研] 材料专硕326求调剂 +6 墨煜姒莘 2026-03-15 7/350 2026-03-17 17:10 by ruiyingmiao
[考研] 考研调剂 +3 淇ya_~ 2026-03-17 5/250 2026-03-17 09:25 by Winj1e
[考研] 070303 总分349求调剂 +3 LJY9966 2026-03-15 5/250 2026-03-16 14:24 by xwxstudy
[考研] 070305求调剂 +3 mlpqaz03 2026-03-14 4/200 2026-03-15 11:04 by peike
[考研] 本科南京大学一志愿川大药学327 +3 麦田耕者 2026-03-14 3/150 2026-03-14 20:04 by 外星文明
信息提示
关闭
请填处理意见
关闭