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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】30金币求一程序
111/2返回列表12下一页
查看: 1065  |  回复: 10
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

njhx505

木虫 (正式写手)

[交流] 【求助】30金币求一程序已有4人参与


请各位高手帮个忙，CA,CB是不同时间测得浓度值，求k1,k2

数据给出，希望高手可以帮忙，谢谢

[ Last edited by njhx505 on 2010-6-4 at 19:42 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2010-06-04 16:22:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lijinfeng042

木虫 (小有名气)

Matlab
★
robert2020(金币+1):鼓励应助！ 2010-06-04 22:30:45
引用回帖:
Originally posted by njhx505 at 2010-06-04 16:22:41:

请各位高手帮个忙，CA,CB是不同时间测得浓度值，求k1,k2

想明确一点
只有两个值?拟合 应该不可以 还是考虑数值积分吧

[ Last edited by lijinfeng042 on 2010-6-4 at 17:19 ]
赞一下(1人) 回复此楼
工作了,偶尔会上来~可以关注新浪微博 @云是风的梦_Matlab
2楼2010-06-04 17:03:55
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师
★
robert2020(金币+1):鼓励应助！ 2010-06-04 22:31:02
上数据

看微分方程式，dCA/dt=-2*dCB/dt

那么只需要知道CA0，CB0 以及 A或B随时间的变化数据就可以估值了

也就是说只需要估值一个微分方程就可以了

[ Last edited by change0618 on 2010-6-4 at 20:12 ]
赞一下(1人) 回复此楼
3楼2010-06-04 17:31:59
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师
zzuwangshilei:积极参与 2010-06-05 10:53:08
数据有的恶劣
回复此楼
4楼2010-06-04 20:28:04
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

njhx505

木虫 (正式写手)
希望高手们能把MATLAB运行程序给我，谢谢
赞一下 回复此楼
5楼2010-06-04 20:39:35
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wudi_82

★ ★
zzuwangshilei(金币+2):辛苦了，多谢应助 2010-06-05 10:52:39
njhx505(金币+3):谢谢 2010-06-06 21:57:51
(2)式*2+(1)式得到   dCA(t)/dt+2*dCB(t)/dt=0;
在任意小区间[t0,t]积分得到：CA(t)+2*CB(t)=const ....(*)式
计算一下你表格中的数据，不满足是常数。至于怎么处理，或者重新采集数据或者像你提到的满足某种最优。
另外：
将(*)式代入到上面其中一个方程比如(1)得到：
dCA/dt=-k1*CA^2-k2*CA+const; 这个方程建议还是数值求解比较方便，比如向前欧拉，向后欧拉方法等等，要想得到高精度的就用Runge-Kutta。
赞一下(1人) 回复此楼
6楼2010-06-05 03:24:37
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

njhx505

木虫 (正式写手)
不能编写一个程序来求解吗，这只是一小部分数据，要是直接手算计算量太大
赞一下 回复此楼
7楼2010-06-05 12:50:54
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师
★
zzuwangshilei(金币+1):多谢积极参与 2010-06-05 15:16:11
njhx505(金币+3): 2010-06-06 21:58:24
你给的CA数据有问题

按照你的微分方程式可以得到一下关系：

CA-CA0=-2*(CB-CB0)

所以我认为你的数据拟合只需要

t=0时的 CA0  CB0
以及一组CA值或者CB值就可以了。

当然CA，CB数据都有的话，可以验证你的数据的准确性

[ Last edited by change0618 on 2010-6-5 at 15:02 ]
赞一下(1人) 回复此楼
8楼2010-06-05 14:58:06
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

njhx505

木虫 (正式写手)
CA0也是知道的，第一排就是CA0,CB0，希望能够帮忙编写个求解k1.k2的程序
赞一下 回复此楼
9楼2010-06-05 16:38:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师
★ ★
kuhailangyu(金币+2):欢迎积极参与 2010-06-05 22:36:53
njhx505(金币+24):谢谢你的程序，真是救了我的急 2010-06-06 21:59:25
引用回帖:
Originally posted by njhx505 at 2010-06-05 16:38:43:
CA0也是知道的，第一排就是CA0,CB0，希望能够帮忙编写个求解k1.k2的程序

给你一个程序，但是数据不是你的，自己修改去吧


function Kinetics
% 动力学ODE方程模型的参数估计
clear all
clc
ExpData = ...
[        0          0.1883
    0.0100    0.2047
    0.0200    0.2181
    0.0300    0.2291
    0.0400    0.2382
    0.0500    0.2459
    0.0600    0.2523
    0.0700    0.2576
    0.0800    0.2622
    0.0900    0.2660
    0.1000    0.2692
    0.1100    0.2719
    0.1200    0.2742
    0.1300    0.2761
    0.1400    0.2777
    0.1500    0.2790
    0.1600    0.2801
    0.1700    0.2811
    0.1800    0.2819
    0.1900    0.2825
    0.2000    0.2830
];
t = ExpData(:,1);   % ExpData第一列为时间
CB = ExpData(:,2);  % ExpData第二列为组分Ｂ的浓度
CB0 = CB(1);        % t=0时，组分B初始浓度
CA0 = 0.8;          % t=0时，组分A初始浓度
CA = CA0-2*(ExpData(:,2)-CB0);  % 由微分方程式导出来的CA与CB关系式
k0 = [20 50];       % 估值参数的猜想值
lb = [0 0];         % 设定的估值参数上限
ub = [+inf  +inf];  % 设定的估值参数下限
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],t,[CA0,CB0],[CA,CB]);
ci = nlparci(k,residual,jacobian);  % 计算置信区间
fprintf('\n\n估计参数值为:\n')
fprintf('\tk1 = %.6f\t置信区间：[%.6f %.6f]\n',k(1),ci(1,: ))
fprintf('\tk2 = %.6f\t置信区间：[%.6f %.6f]\n',k(2),ci(2,: ))
tt = linspace(t(1),t(end),101);
[tt C] = ode45(@KineticEqs,tt,[CA0,CB0],[],k);  % 由估值k计算时间序列tt下的A，B浓度
figure(1)
plot(t,CA,'o',tt,C(:,1),'r-')
xlabel('t'); ylabel('C_A')
figure(2)
plot(t,CB,'o',tt,C(:,2),'r-')
xlabel('t'); ylabel('C_B')
% ------------------------------------------------------------------
function f = ObjFunc(k,tspan,x0,yexp)
[t,y] = ode45(@KineticEqs,tspan,x0,[],k);
f = y(: ) - yexp(: );
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
dxdt =[-k(1)*x(1)^2+2*k(2)*x(2);
       0.5*k(1)*x(1)^2-k(2)*x(2)];

[ Last edited by change0618 on 2010-6-5 at 18:23 ]
赞一下(1人) 回复此楼
10楼2010-06-05 18:15:24
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 njhx505 的主题更新
111/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭