24小时热门版块排行榜    

查看: 3388  |  回复: 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的回帖

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的回帖
查看全部 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的回帖

皓小天

木虫之王 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
我记得已经回复了你这个问题,http://muchong.com/bbs/viewthread.php?tid=9255377&authorid=1864745
请不要站内找我要书,如果需要请到书籍板块求助
4楼2015-08-31 22:03:03
已阅   回复此楼   关注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的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 我的国基提前知道中了,可是同事的操作让我实在接受不了,怎么会有这样的人 +9 家与远方 2026-08-10 14/700 2026-08-13 09:41 by 臭臭不臭01
[基金申请] 结合人工智能,周易传统文化,filecode打分制来了,3分以上希望很大。 +3 Tide man 2026-08-12 4/200 2026-08-13 08:35 by ZJTJZ
[基金申请] 应该是93bebmhtak前后十一个字符比较关键 +22 Lanmanbaby 2026-08-09 36/1800 2026-08-12 22:50 by sdfapple719
[基金申请] 好奇怪的filecode +5 布布和一二 2026-08-08 6/300 2026-08-12 16:11 by 云上清扬
[基金申请] filecode +13 documentary 2026-08-10 15/750 2026-08-12 15:56 by 云上清扬
[论文投稿] 职称评审,求友友推荐见刊最快的期刊 +4 工厂打螺丝 2026-08-08 4/200 2026-08-12 09:18 by hansi2025
[硕博家园] 读博的好处 +3 lnee 2026-08-11 3/150 2026-08-11 18:10 by 希望我好好的
[基金申请] 什么时候出结果,有咨询渠道??? +3 Tide man 2026-08-11 3/150 2026-08-11 17:54 by kudofaye
[基金申请] 听说今天filecode变了 +26 布布和一二 2026-08-06 49/2450 2026-08-11 13:18 by WH3796
[基金申请] 基金中了 +15 laoda193707 2026-08-06 15/750 2026-08-11 00:11 by jiafei2190
[基金申请] 静等基金结果 +5 gjjjzhong 2026-08-10 16/800 2026-08-10 17:19 by Tide man
[基金申请] 据悉今年马上要出结果了 +7 瞬息宇宙 2026-08-10 8/400 2026-08-10 12:42 by Vivilian
[基金申请] 2026国自然放榜时间 +9 布布和一二 2026-08-08 9/450 2026-08-10 11:22 by xxxx2020
[基金申请] 这样的filecode谁见过 +11 布布和一二 2026-08-08 22/1100 2026-08-10 11:10 by wmfsnow
[基金申请] 面上项目filecode邪修 +5 西山十月 2026-08-09 7/350 2026-08-10 07:32 by 仁砚薪传
[基金申请] filecode与中标关系的预测 +5 布布和一二 2026-08-07 5/250 2026-08-09 16:15 by 袁向阳007
[基金申请] 关于filecode,很负责任的告诉大家 +6 爱看书的可乐 2026-08-08 7/350 2026-08-08 22:13 by a_niu
[基金申请] 关于豆爷回答的JTJC与%2F数量 +5 yang182083 2026-08-06 7/350 2026-08-08 18:29 by zhanghaozhu
[基金申请] 化学口download_prp&fileCode的固定段好像这几天一直没变,有变的大神么? +3 Tide man 2026-08-07 4/200 2026-08-07 22:39 by Tide man
[基金申请] filecode +14 等待解的谜 2026-08-06 19/950 2026-08-07 12:20 by wlwhappy
信息提示
请填处理意见