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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 微分方程ode45求解，最小二乘法优化微分方程参数，程序运行求助
61/1返回列表
查看: 2118  |  回复: 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的回帖

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的回帖

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的回帖

whqs8426212

铜虫 (正式写手)

【答案】应助回帖

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

[ 发自小木虫客户端 ]
赞一下 回复此楼
5楼2014-09-20 18:18:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

birdcanfly

金虫 (正式写手)
引用回帖:
5楼: Originally posted by whqs8426212 at 2014-09-20 18:18:39
两个边界值范围缩小，试试

如果不知道参数值，如何缩小边界？刚才试了试退火算法，也是一样模拟值和给出的初始值关系密切。
赞一下 回复此楼
6楼2014-09-23 18:08:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 birdcanfly 的主题更新
61/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 放榜前的不淡定 40+4 snowwithsea 2026-08-19 13/650 2026-08-21 08:30 by 北京莱茵编辑
[基金申请] 让我中一个面上吧！ +9 大萍1987 2026-08-20 11/550 2026-08-21 08:19 by littlezl
[基金申请] 估计是周四 +3 archvillain 2026-08-18 3/150 2026-08-21 01:48 by jnhyjjm
[基金申请] 基金啊基金 +3 longfie172 2026-08-20 3/150 2026-08-21 00:04 by wuchongjun
[基金申请] 科研孤儿太难了 +14 我4大白菜 2026-08-20 15/750 2026-08-20 22:33 by 袁向阳007
[基金申请] 言之凿凿，明天官网会通知放榜？ +15 医学老男孩 2026-08-20 17/850 2026-08-20 20:35 by CaoLG
[基金申请] 今天系统多次维护，明天很可能放榜！ +10 zju2000 2026-08-16 11/550 2026-08-20 20:02 by cl479861084
[基金申请] 时间戳变了，能看出什么问题？ +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[基金申请] filecode，4个jtjc了 +11 ziyangfang 2026-08-19 13/650 2026-08-19 21:28 by aasahr
[教师之家] 为什么余额宝的年化利率越来越低？主要原因有哪些？ +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] filecode=后面第一个是大写字母 +10 wangze12014 2026-08-14 12/600 2026-08-19 16:56 by 苦难博士
[基金申请] 2027广东省杰青 +4 奶牛小黑 2026-08-15 10/500 2026-08-19 11:02 by wanfengnew
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 快农历七夕节了，轻松一下，男人悄悄话，女施主请不要进来。 +6 Tide man 2026-08-14 7/350 2026-08-19 09:56 by ZJTJZ
[基金申请] 什么时候开奖？ +6 CrisMessi 2026-08-18 6/300 2026-08-19 08:06 by Equinoxhua
[基金申请] 明天放榜？ +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[论文投稿] 投稿咨询 +4 wwm09 2026-08-17 6/300 2026-08-18 15:36 by wwm09
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[精细化工] 招聘 金属平磨液，抛光液研发工程师 +3 小天0311 2026-08-14 3/150 2026-08-16 07:31 by H9PLUS
[基金申请] 各位道友，我要去昆明玩几天，回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
信息提示
关闭
请填处理意见
关闭