24小时热门版块排行榜    

查看: 3253  |  回复: 5
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

zhl717

铜虫 (小有名气)

[求助] matlab拟合微分方程组中的参数 已有2人参与

最近一直在看matlab求解微分方程组中参数问题,但是由于本人以前没有接触过matlab,所以还是不太懂,请懂的大神帮忙写一下,谢谢!!

有高版本1st opt的大神也可以帮忙跑一下,上回发了一个帖子,没人回应,谢谢各位!!

dx/dt=a*x*(1-x/b)*(1+s/c)^(-1)
dp/dt=d*dx/dt+e*[s/(s+f)]*x
ds/dt=g*dx/dt+h*dp/dt+i*x

t为自变量,x、p、s为因变量,a-i  均为参数

数据如下:
t,   X,   P,       S
0   0.35  0      99.8
12  2.75  0.55   88.96
24  4.6   3.72   73.49
36  6.27  8.19   61.57
48  8.07  14.12  47.38
60  9.73 19.37  33.63
72  10.41 22.54  24.65
84  10.7  27.61  11.19
96  10.53 32.49  3.2
108 10.56 35.62  0
120 10.59 37.58  0
132 10.4  38.44  0
144 10.77 39.88  0
回复此楼

» 本帖@通知

» 猜你喜欢

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

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

narian

木虫 (著名写手)

总司令

引用回帖:
2楼: Originally posted by 月只蓝 at 2015-08-31 16:48:55
要拟合九个参数,这种问题还是推荐用1stOpt,参数个数太多,用MATLAB比较困难。
下面给出MATLAB参考代码,代码中初值k0可以试着调一下,我初步调了一下初值,拟合效果不好。

function k1k2k3
format long
cle ...

高手,有问题请教

[ 发自小木虫客户端 ]
回不去的青春,那是我一生一次的认真!
5楼2015-09-01 06:30:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +1
zhl717: 金币+5, ★★★很有帮助 2015-08-31 17:03:01
要拟合九个参数,这种问题还是推荐用1stOpt,参数个数太多,用MATLAB比较困难。
下面给出MATLAB参考代码,代码中初值k0可以试着调一下,我初步调了一下初值,拟合效果不好。
CODE:
function k1k2k3
format long
clear all
clc
tspan = 0:12:144;
x0 = [0.35  0      99.8];

k0 = [0.5 8 1.5 2 -0.5 0.1 0.5 5 -0.01];   %k的初值,最需要调节

lb = [1 1 1 1 1 1 1 1 1 ]*-inf;
ub = [1 1 1 1 1 1 1 1 1]*inf;

data=...
    [


12  2.75  0.55   88.96
24  4.6   3.72   73.49
36  6.27  8.19   61.57
48  8.07  14.12  47.38
60  9.73 19.37  33.63
72  10.41 22.54  24.65
84  10.7  27.61  11.19
96  10.53 32.49  3.2
108 10.56 35.62  0
120 10.59 37.58  0
132 10.4  38.44  0
144 10.77 39.88  0
];
yexp = data(:,2:4);
options=optimset('MaxFunEvals',1500);

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);

[t uu]=ode45(@KineticsEqs,tspan,x0,[],k);
uexp=[0   0.35  0      99.8
12  2.75  0.55   88.96
24  4.6   3.72   73.49
36  6.27  8.19   61.57
48  8.07  14.12  47.38
60  9.73 19.37  33.63
72  10.41 22.54  24.65
84  10.7  27.61  11.19
96  10.53 32.49  3.2
108 10.56 35.62  0
120 10.59 37.58  0
132 10.4  38.44  0
144 10.77 39.88  0
];
figure(1)

size(uexp(:,2));
size(uexp(:,3));
size(uexp(:,4));
size(uu(:,1));
size(uu(:,2));
size(uu(:,3));

R1_Square=1-sum((uexp(2:end,2)-uu(2:end,1)).^2)./sum((uexp(2:end,2)-mean( uu(2:end,1)  )).^2);
R2_Square=1-sum((uexp(2:end,3)-uu(2:end,2)).^2)./sum((uexp(2:end,3)-mean( uu(2:end,2)  )).^2);
R3_Square=1-sum((uexp(2:end,4)-uu(2:end,3)).^2)./sum((uexp(2:end,4)-mean( uu(2:end,3)  )).^2);



plot(tspan,uexp(:,2),'or',tspan,uu(:,1),'-r',tspan,uexp(:,3),'>b',tspan,uu(:,2),'b-',tspan,uexp(:,4),'*k',tspan,uu(:,3),'k-')


fprintf('\n\tk1=%.9f',k(1))
fprintf('\n\tk2=%.9f',k(2))
fprintf('\n\tk3=%.9f',k(3))
fprintf('\n\tk4=%.9f',k(4))
fprintf('\n\tk5=%.9f',k(5))
fprintf('\n\tk6=%.9f',k(6))
fprintf('\n\tk7=%.9f',k(7))
fprintf('\n\tk8=%.9f',k(8))
fprintf('\n\tk9=%.9f',k(9))

fprintf('\n\t方程一的决定系数R1_Square=%.6f',R1_Square);
fprintf('\n\t方程二的决定系数R2_Square=%.6f',R2_Square);
fprintf('\n\t方程三的决定系数R3_Square=%.6f',R3_Square);




function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,2);

ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);


f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3))];

function dCdt = KineticsEqs(t,u,k)  
x=u(1);
p=u(2);
s=u(3);

a=k(1);
b=k(2);
c=k(3);
d=k(4);
e=k(5);
f=k(6);
g=k(7);
h=k(8);
i=k(9);



dxdt = a*x*(1-x/b).*(1+s/c).^(-1);
dpdt=d*dxdt+e*(s./(s+f)).*x;
dsdt=g*dxdt+h*dpdt+i*x;
dCdt = [dxdt;dpdt;dsdt];

计算结果:
        k1=-24.867032689
        k2=5.968347685
        k3=0.525742903
        k4=11.267708389
        k5=6.039602224
        k6=-1.614606972
        k7=-15.309885669
        k8=6.596877645
        k9=0.664727164
        方程一的决定系数R1_Square=-0.002145
        方程二的决定系数R2_Square=0.359295
        方程三的决定系数R3_Square=-0.035619
matlab拟合微分方程组中的参数
附图1.png

MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
2楼2015-08-31 16:48:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhl717

铜虫 (小有名气)

引用回帖:
2楼: Originally posted by 月只蓝 at 2015-08-31 16:48:55
要拟合九个参数,这种问题还是推荐用1stOpt,参数个数太多,用MATLAB比较困难。
下面给出MATLAB参考代码,代码中初值k0可以试着调一下,我初步调了一下初值,拟合效果不好。

function k1k2k3
format long
cle ...

好的,非常感谢
加油!!!
3楼2015-08-31 17:03:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

皓小天

木虫之王 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
我记得已经回复了你这个问题,http://muchong.com/bbs/viewthread.php?tid=9255377&authorid=1864745
请不要站内找我要书,如果需要请到书籍板块求助
4楼2015-08-31 22:03:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 305分求调剂(食品工程) +5 Sxy112 2026-03-21 7/350 2026-03-24 12:27 by 544594351
[考研] 一志愿哈工大,085400,320,求调剂 +3 gdlf9999 2026-03-24 3/150 2026-03-24 10:08 by 搏击518
[考研] 一志愿河北工业大学0817化工278分求调剂 +7 jhybd 2026-03-23 12/600 2026-03-24 09:03 by jhybd
[基金申请] 请教下大家 2026年国家基金申请是双盲审吗? +3 lishucheng1 2026-03-22 5/250 2026-03-24 08:22 by gltch
[考研] 材料专硕找调剂 +3 哈哈哈吼吼吼哈 2026-03-23 3/150 2026-03-23 23:13 by peike
[考研] 一志愿陕师大生物学071000,298分,求调剂 +3 SYA! 2026-03-23 3/150 2026-03-23 19:09 by macy2011
[考研] 生物学一志愿985,分数349求调剂 +6 zxts12 2026-03-21 9/450 2026-03-23 18:37 by macy2011
[考研] 276求调剂 +3 YNRYG 2026-03-21 4/200 2026-03-23 08:31 by 醉在风里
[考研] 0854电子信息求调剂 +3 α____ 2026-03-22 3/150 2026-03-22 21:28 by zhq0425
[考研] 354求调剂 +7 Tyoumou 2026-03-18 10/500 2026-03-22 11:11 by 人来盛
[考博] 招收博士1-2人 +3 QGZDSYS 2026-03-18 4/200 2026-03-22 10:25 by QGZDSYS
[考研] 初试 317 +7 半拉月丙 2026-03-20 7/350 2026-03-21 22:26 by peike
[考研] 求调剂 +4 要好好无聊 2026-03-21 4/200 2026-03-21 18:57 by 学员8dgXkO
[考研] 一志愿南大,0703化学,分数336,求调剂 +3 收到VS 2026-03-21 3/150 2026-03-21 18:42 by 学员8dgXkO
[考研] 297求调剂 +11 戏精丹丹丹 2026-03-17 12/600 2026-03-21 17:47 by ColorlessPI
[考研] 化学求调剂 +4 临泽境llllll 2026-03-17 5/250 2026-03-21 02:23 by JourneyLucky
[考研] 085700资源与环境308求调剂 +12 墨墨漠 2026-03-18 13/650 2026-03-21 01:42 by JourneyLucky
[考研] 考研调剂求学校推荐 +3 伯乐29 2026-03-18 5/250 2026-03-20 22:59 by JourneyLucky
[考研] 0856调剂,是学校就去 +8 sllhht 2026-03-19 9/450 2026-03-20 14:25 by 无懈可击111
[论文投稿] 申请回稿延期一个月,编辑同意了。但系统上的时间没变,给编辑又写邮件了,没回复 10+3 wangf9518 2026-03-17 4/200 2026-03-19 23:55 by babero
信息提示
请填处理意见