查看: 6559  |  回复: 21

yupeifeng

金虫 (小有名气)

[交流] matlab求解方程中的参数已有4人参与

搞了这么久还是没有解决这个问题。现在求助于高手给予指点。谢谢了!!
实验数据为:(t,c)=(0,0.69)(2,0.645)(4,0.635)(8,0.62)(24,0.61)(48,0.61).其中t为时间,c为某离子的浓度。
动力学方程模型为:-dc/dt=k*(c0-c)^(1/3)*(c-c~).
其中c0为初始浓度可以取0.7,c~为平衡浓度取0.61.
怎么样才能拟合出参数k的值呢?谢谢大家给出程序代码,再次感谢
回复此楼

» 猜你喜欢

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

憧憬未来
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

dbb627

荣誉版主 (著名写手)

★ ★ ★ ★ ★
小木虫: 金币+0.5, 给个红包,谢谢回帖
fegg7502: 金币+4, 多谢交流 2012-07-09 08:07:15
这个其实论坛里有很多例子,参考就能写出来。
给你写了个
CODE:
function parafit
clear all;
t=[0 2 4 8 24 48];
y=[0.69 0.645 0.635 0.62 0.61 0.61];
y0=0.69;

% Nonlinear least square estimate using lsqnonlin()
beta0=0.5;
lb=[0];ub=[inf];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@Func,beta0,lb,ub,[],t,y,y0);         
ci = nlparci(beta,residual,jacobian);
beta;
% result
fprintf('\n Estimated Parameters by Lsqnonlin():\n')
fprintf('\t k = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('  The sum of the residual squares is: %.1e\n\n',sum(residual.^2))

% plot of fit results
tspan = [0  max(t)];
[tt yc] = ode45(@ModelEqs,tspan,y0,[],beta);
tc=linspace(0,max(t),200);
yca = spline(tt,yc,tc);
plot(t,y,'ro',tc,yca,'r-');
hold on
xlabel('Time');
ylabel('Concentration');
hold off
% =======================================
function f1 = Func(beta,t,y,y0)        % Define objective function
tspan =t;
[tt yy] = ode45(@ModelEqs,tspan,y0,[],beta);
yc= spline(tt,yy,t);
f1=y-yc;
% ==================================
function dydt = ModelEqs(t,y,beta)          % Model equations
dydt = -beta*(0.7-y).^(1/3)*(y-0.61);

The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
2楼2012-07-02 17:38:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

dbb627

荣誉版主 (著名写手)

★ ★ ★ ★
fegg7502: 金币+4, 应助指数+1, 专家考核, 多谢交流 2012-07-09 08:07:38
结果如下
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to
its initial value is less than the default value of the function tolerance.




Estimated Parameters by Lsqnonlin():
         k = 1.0138 ± 0.3000
  The sum of the residual squares is: 1.5e-004


The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
3楼2012-07-02 17:39:27
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yupeifeng

金虫 (小有名气)

引用回帖:
3楼: Originally posted by dbb627 at 2012-07-02 17:39:27
结果如下
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to
its initial value is less than the default value of the function tolerance.

< ...

首先感谢这位 学长啊,我有点看不懂啊,你帮我 分析下思路,步骤是什么?就是先做设么计算,
憧憬未来
4楼2012-07-02 22:01:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

★ ★ ★
小木虫: 金币+0.5, 给个红包,谢谢回帖
fegg7502: 金币+2, 多谢交流 2012-07-09 08:07:55
1stOpt求解:
CODE:
Constant c0=0.7,cp=0.61;
Variable t,c;
ODEFunction c'=-k*(c0-c)^(1/3)*(c-cp)
Data;
t,c
0,0.69
2,0.645
4,0.635
8,0.62
24,0.61
48,0.61

均方差(RMSE): 0.00535058522113925
残差平方和(SSE): 0.000143143811043369
相关系数(R): 0.955733024504782
相关系数之平方(R^2): 0.913425614129058
决定系数(DC): 0.852429060780033

参数                  最佳估算
--------------------        -------------
k        1.02159845408008
5楼2012-07-03 08:52:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yupeifeng

金虫 (小有名气)

引用回帖:
5楼: Originally posted by dingd at 2012-07-03 08:52:49
1stOpt求解:

Constant c0=0.7,cp=0.61;
Variable t,c;
ODEFunction c'=-k*(c0-c)^(1/3)*(c-cp)
Data;
t,c
0,0.69
2,0.645
4,0.635
8,0.62
24,0.61
48,0.61


均方差(RMSE): 0.00535058522113925
...

你好,我今天下载了个1stopt,不怎么会用啊,这个和matlab一样吗?
憧憬未来
6楼2012-07-03 09:31:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

★ ★
小木虫: 金币+0.5, 给个红包,谢谢回帖
fegg7502: 金币+1, 多谢交流 2012-07-09 08:08:09
比Matlab简单啊,优化拟合方面更是远胜Matlab!
7楼2012-07-03 09:45:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yupeifeng

金虫 (小有名气)

引用回帖:
7楼: Originally posted by dingd at 2012-07-03 09:45:41
比Matlab简单啊,优化拟合方面更是远胜Matlab!

不太会用,怎么去定义哪些量,怎么运行呢?
憧憬未来
8楼2012-07-03 09:49:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yupeifeng

金虫 (小有名气)

引用回帖:
5楼: Originally posted by dingd at 2012-07-03 08:52:49
1stOpt求解:

Constant c0=0.7,cp=0.61;
Variable t,c;
ODEFunction c'=-k*(c0-c)^(1/3)*(c-cp)
Data;
t,c
0,0.69
2,0.645
4,0.635
8,0.62
24,0.61
48,0.61


均方差(RMSE): 0.00535058522113925
...

我把你的程序复制后,提示没有定义参数,这是怎么回事啊
憧憬未来
9楼2012-07-03 10:16:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

★ ★
小木虫: 金币+0.5, 给个红包,谢谢回帖
fegg7502: 金币+1, 多谢交流 2012-07-09 08:08:30
你用的是网上的1.5PJ版吧?老了点,好像也不支持微分方程计算。用新版本没问题。
10楼2012-07-03 10:38:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 yupeifeng 的主题更新
普通表情 高级回复(可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[育儿交流] 分享一个蒜皮小事 +17 lxy12365492 2024-11-17 45/2250 2024-11-24 00:36 by gxt273
[硕博家园] 换导师学院里面不同意咋办,原导师现导师都同意了 +12 成成的老婆 2024-11-22 18/900 2024-11-23 23:28 by 单纯弱小
[论文投稿] 投稿返修中,增加了作者,论文被拒了,还有办法补救嘛? +9 mollyzhang_2003 2024-11-22 18/900 2024-11-23 22:02 by mollyzhang_2003
[考博] 求个博导把我捞走,一篇JMCA,另有一篇AFM共一 +5 骑萝卜钓鱼 2024-11-21 7/350 2024-11-23 21:03 by Zixiao_Wang
[硕博家园] 搞科研是不是很锻炼身体?高血压,痛风是小病 +4 akslis2024 2024-11-23 8/400 2024-11-23 20:51 by 328838485
[基金申请] 长江学者评审开始了,有评审专家爆料被地毯式打招呼 +13 babu2015 2024-11-22 14/700 2024-11-23 14:32 by ca0yan9
[考博] 读博求助 10+3 姜大胖仙女 2024-11-18 5/250 2024-11-23 03:53 by dosupy
[教师之家] 坐标西安,找个搭子 +10 hamihappy 2024-11-19 12/600 2024-11-22 23:14 by yiran909
[考博] 招申请考核博士研究生 +4 zhangqi1012 2024-11-22 4/200 2024-11-22 22:24 by 求上岸的
[基金申请] 美国德州州立高校将禁止和中国合作科研 +5 babu2015 2024-11-20 5/250 2024-11-22 21:16 by jurkat.1640
[基金申请] 2024 年国社科后期资助项目到什么阶段了? +6 strong409 2024-11-17 8/400 2024-11-22 18:53 by bnugjsyh
[考博] 24级双9硕钙钛矿材料方向求博导收留 +6 yfdsefx 2024-11-19 8/400 2024-11-22 14:54 by 鱼儿梦水深
[论文投稿] 审稿人让补实验可以拒绝吗? 28+6 maxlhy 2024-11-20 13/650 2024-11-22 12:07 by suifengdao
[有机交流] 实验室减压蒸馏二甲苯 20+4 lijiajun09 2024-11-19 8/400 2024-11-21 22:50 by s66230911
[论文投稿] assigning reviewer一周就awaiting ECI decision是什么意思 +3 mollyzhang_2003 2024-11-21 5/250 2024-11-21 20:08 by mollyzhang_2003
[论文投稿] ACS Catalysis投稿问题求助解答 10+5 葛旺鑫 2024-11-18 5/250 2024-11-21 14:43 by TopEdit
[论文投稿] under review 50+7 星辰2022 2024-11-19 13/650 2024-11-21 14:39 by TopEdit
[论文投稿] 返修意见:有关于XRD精修R值过高的问题 5+4 zychiang 2024-11-18 10/500 2024-11-20 19:33 by zychiang
[论文投稿] 请问JCR分区与WOS分区是一样的吗? 5+4 zhongyuan764 2024-11-19 7/350 2024-11-20 14:19 by cqupenghao
[有机交流] NBS溴代 100+3 风月客 2024-11-17 9/450 2024-11-19 18:15 by 风月客
信息提示
请填处理意见