版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

雾隐村的白

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 2155.8
帖子: 30
在线: 17.1小时
虫号: 3688807
注册: 2015-02-14
专业: 高分子合成化学

[求助] matlab模拟微分方程参数求助各位大佬已有1人参与

反应动力学方程：dC/dt=-k1*c*(2.413+C)+k2*(4.826-C)^2;
t=[0 15 30 45 60 75 90 120];
C=[4.826 4.206045728 3.081681077 2.582976758 2.368099268 2.296997119 2.259547446 2.221752483];
求解参数 k1,k2；
1stopt没有软件条件，所以只能寄希望于matlab。
请教各位大佬matlab代码怎么写，感激不尽！！！

回复此楼

» 猜你喜欢

反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有65人回复
新西兰Robinson研究所招收全奖PhD 已经有0人回复
石墨烯转移--二氧化硅衬底石墨烯已经有0人回复
笼目材料中量子自旋液体基态的证据已经有0人回复
数学教学论硕士可以读数学物理博士吗？已经有0人回复
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生已经有4人回复
澳门大学应用物理及材料工程研究院潘晖教授课题组诚招博士后已经有11人回复

1楼 2019-12-16 09:15:36

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

独孤神宇

版主 (知名作家)

应助: 490 (硕士)
贵宾: 0.008
金币: 31014.8
散金: 802
红花: 122
沙发: 1
帖子: 5600
在线: 855.5小时
虫号: 3522474
注册: 2014-11-06
性别: GG
专业: 机械动力学
管辖: 计算模拟

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
感谢参与，应助指数 +1
雾隐村的白: 金币+100, ★★★很有帮助 2019-12-16 13:21:59

http://blog.sina.com.cn/s/blog_c0cb8ce60102ysqt.html

赞一下(1人)

回复此楼

» 本帖已获得的红花（最新10朵）

雾隐村的白

数值计算

2楼2019-12-16 10:21:29

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

雾隐村的白

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 2155.8
帖子: 30
在线: 17.1小时
虫号: 3688807
注册: 2015-02-14
专业: 高分子合成化学

送红花一朵

引用回帖:

2楼: Originally posted by 独孤神宇 at 2019-12-16 10:21:29
http://blog.sina.com.cn/s/blog_c0cb8ce60102ysqt.html

模仿代码修改如下：
function ODEfunction
clear all;clc
format long
tspan=[0 15 30 45 60 75 90 120]; % size=1*8
yexp=[4.826 4.206045728 3.081681077 2.582976758 2.368099268 2.296997119 2.259547446 2.221752483]'; % size=8*1,已将第一个数据取出作为下面的初始值
k0=[0.002 0.003]; %猜测初值
y0=4.826;          % 初始状态
lb=[0 0 ];          % 参数下限
ub=[1 1];    % 参数上限
yy=[y0 yexp'];       % 未给定初始值，则第一行作为初始值
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],y0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;

% 这一步通常被省略，通过反复迭代初始值得到最优解，加上后可以降低对初始值的依赖。
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
% 需要指出，这种方法并非在所有场合均有效，但有时确实可以改善求解效果。

k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],y0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('The sum of the squares is: %.1e\n\n',resnorm)
%---------------------------------------------------------------------
ts=0:0.5:max(tspan);       %用于计算的步长，步数可比实际数据多
[ts,ys]=ode45(@KineticEqs,ts,y0,[],k);       %微分方程求解
[ttt,XXsim] = ode45(@KineticEqs,tspan,y0,[],k); %指定点微分方程求解
y=XXsim(2:end);                   % 与实际数数据维数保持一致
R2=1-sum((yexp-y).^2)./sum((yexp-mean(y)).^2);
fprintf('\n\t决定系数R-Square = %.6f',R2);
figure
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best');
xlabel('时间');ylabel('计算结果');
% ------------------------------------------------------------------
function f = ObjFunc4Fmincon(k,x0,yexp)
tspan = 0 : 1 : 8;                         % ts=0:1:max(tspan);
[t,Xsim] = ode45(@KineticEqs,tspan,x0,[],k);  % ode45函数参数传递的调用形式
y = Xsim(2:end);                            % 对应实验数据  yexp
f = sum((y-yexp).^2);                      % 计算平方和，供fmincon调用
%---------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)          % lsqnonlin目标函数
tspan = 0: 1 : 8;
[t,Xsim] = ode45(@KineticEqs,tspan,x0,[],k);
ysim = Xsim(2:end);                         % 确保维数一致
f=ysim-yexp;
%----------------------------------------------------------
function dydt = KineticEqs(t,y,k)          % 微分方程
beta(1)=k(1);
beta(2)=k(2);
dydt = -beta(1)*y*(2.413+y)+beta(2)*(4.826-y)^2;
%code end

---------------------------------------------------------------------
使用函数fmincon()估计得到的参数值为:
k1 = 0.0203
k2 = 0.0021
The sum of the squares is: 9.6e-01
出现了一些错误：
Error using  -
Matrix dimensions must agree.

Error in ODEfunction (line 36)
R2=1-sum((yexp-y).^2)./sum((yexp-mean(y)).^2);
也没有拟合的图表出现。

您有时间看看是什么问题么？

赞一下

回复此楼

3楼2019-12-16 11:28:42

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

独孤神宇

版主 (知名作家)

应助: 490 (硕士)
贵宾: 0.008
金币: 31014.8
散金: 802
红花: 122
沙发: 1
帖子: 5600
在线: 855.5小时
虫号: 3522474
注册: 2014-11-06
性别: GG
专业: 机械动力学
管辖: 计算模拟

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
雾隐村的白: 金币+100, ★★★★★最佳答案 2019-12-16 14:34:45

引用回帖:

3楼: Originally posted by 雾隐村的白 at 2019-12-16 11:28:42
模仿代码修改如下：
function ODEfunction
clear all;clc
format long
tspan=; % size=1*8
yexp='; % size=8*1,已将第一个数据取出作为下面的初始值
k0=; %猜测初值
y0=4.826; % 初始状态
...

function ODEfunction_12_16
clear all;clc
format long
tspan=[0 15 30 45 60 75 90 120]; % size=1*8
yexp=[4.826 4.206045728 3.081681077 2.582976758 2.368099268 2.296997119 2.259547446 2.221752483]'; % size=8*1,已将第一个数据取出作为下面的初始值
k0=[0.002 0.003];    %猜测初值
y0=4.826;             % 初始状态
lb=[0 0 ];          % 参数下限
ub=[1 1];             % 参数上限

[k,resnorm] =lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],y0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('The sum of the squares is: %.1e\n\n',resnorm)
%---------------------------------------------------------------------
ts=0:1:max(tspan);                                     %用于计算的步长，步数可比实际数据多
[ts,ys]=ode45(@KineticEqs,ts,y0,[],k);                %微分方程求解
[ttt,XXsim] = ode45(@KineticEqs,tspan,y0,[],k);       %指定点微分方程求解
R2=1-sum((yexp-XXsim).^2)./sum((yexp-mean(XXsim)).^2);
fprintf('\n\t决定系数R^2 = %.6f',R2);
figure
plot(ts,ys,'b',tspan,yexp,'or'),legend('计算值','实验值','Location','best');
xlabel('时间');ylabel('计算结果');
% ------------------------------------------------------------------

%---------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)          % lsqnonlin目标函数
[t,Xsim] = ode45(@KineticEqs,tspan,x0,[],k);
f=Xsim-yexp;
end
%----------------------------------------------------------
function dydt = KineticEqs(t,y,k)          % 微分方程
beta(1)=k(1);
beta(2)=k(2);
dydt = -beta(1)*y*(2.413+y)+beta(2)*(4.826-y)^2;
end
end

赞一下

回复此楼

» 本帖已获得的红花（最新10朵）

雾隐村的白

数值计算

4楼2019-12-16 13:51:12

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

雾隐村的白

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 2155.8
帖子: 30
在线: 17.1小时
虫号: 3688807
注册: 2015-02-14
专业: 高分子合成化学

送红花一朵

引用回帖:

4楼: Originally posted by 独孤神宇 at 2019-12-16 13:51:12
function ODEfunction_12_16
clear all;clc
format long
tspan=; % size=1*8
yexp='; % size=8*1,已将第一个数据取出作为下面的初始值
k0=;    %猜测初值
y0=4.826;             % 初始状态
lb=;    ...

非常感谢！帮了大忙！！！

赞一下

回复此楼

5楼2019-12-16 14:34:31

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