| 查看: 2132 | 回复: 5 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[求助]
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代码怎么写，感激不尽！！！ |
回复此楼» 猜你喜欢
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有122人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有11人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有12人回复
-
桂林理工大学物理学专业招收调剂
已经有13人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
美国顶级物理期刊《应用物理快报》（APL）的编辑欺骗和歧视作者及AIP的官僚主义傲慢
已经有4人回复
送红花一朵 |
模仿代码修改如下： 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
独孤神宇
版主 (知名作家)
- 应助: 490 (硕士)
- 贵宾: 0.008
- 金币: 31016.3
- 散金: 802
- 红花: 122
- 沙发: 1
- 帖子: 5600
- 在线: 856.5小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
感谢参与，应助指数 +1
雾隐村的白: 金币+100, ★★★很有帮助 2019-12-16 13:21:59
感谢参与，应助指数 +1
雾隐村的白: 金币+100, ★★★很有帮助 2019-12-16 13:21:59
» 本帖已获得的红花（最新10朵）
2楼2019-12-16 10:21:29
独孤神宇
版主 (知名作家)
- 应助: 490 (硕士)
- 贵宾: 0.008
- 金币: 31016.3
- 散金: 802
- 红花: 122
- 沙发: 1
- 帖子: 5600
- 在线: 856.5小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
雾隐村的白: 金币+100, ★★★★★最佳答案 2019-12-16 14:34:45
雾隐村的白: 金币+100, ★★★★★最佳答案 2019-12-16 14:34:45
|
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
5楼2019-12-16 14:34:31
