24小时热门版块排行榜    

查看: 1361  |  回复: 0
【悬赏金币】回答本帖问题,作者757272131将赠送您 10 个金币

757272131

新虫 (初入文坛)

[求助] 使用matlab最小二乘法拟合求解常微分方程组未知参数

function k1k2k32
format long
clear all
clc
tspan = [0 6 24 44 68 72 74 92 104 116]';%%这是时间
yexp= [3.111,3.639,3.887,4.289,4.658,5.531,6.218,6.979,7.111,7.114]';%%%这是菌落总数
x0 = [3.111];
k0 = [0.1  3  45 4 37 8.1];
lb = [0  0  0  0 0  0];
ub = [100  100  100 100 100 100];
%
% % opts = statset('nlinfit');
% % opts.RobustWgtFun = 'bisquare';
% % mdll = fitnlm(tspan,yexp,ObjFunc,b0,'Options')
% mdll = fitnlm(tspan,yexp,ObjFunc,b0,'Options',opts)
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);  

ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.9f ± %.9f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.9f ± %.9f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.9f ± %.9f\n',k(6),ci(6,2)-k(6))
fprintf('  The sum of the squares is: %.9e\n\n',resnorm)


tsa=0:0.01:max(tspan);
[tsa ysa]=ode45(@KineticsEqs,tsa,x0,[],k);

figure(1),
plot(tsa,ysa(:,1),'b',tspan,yexp,'or'),legend('计算值','实验值','Location','best');

function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[tspan  y] = ode45(@KineticsEqs,tspan,x0,[],k);
f = (y-yexp);

function dCdt = KineticsEqs(tspan ,C,k)   
T=4*(tspan >=0 & tspan <=30)+10*(tspan >31 & tspan <=90)+4*(tspan >90 ); % 假设这是我写的温度关于时间的分段函数
umax=((k(2)*(T-k(3))*(T-k(4))*(T-k(4)))/(((k(5)-k(4))*(T-k(5))-(k(5)-k(3))*(k(5)+k(4)-2*T))*(k(5)-k(4))));
dCdt=(1/(1+exp(-4*(tspan -k(1)))))*umax*(1-exp(C-k(6)));

这是根据他人的代码更改的,得到的结果是置信区间过大?想有什么办法解决?以及想用fitnlm函数进行拟合,达到一些参数评价t值p值等,但用fitnlm运行不起来。 @beefly
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 757272131 的主题更新
不应助 确定回帖应助 (注意:应助才可能被奖励,但不允许灌水,必须填写15个字符以上)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 QTy3jDtz1uLt 2026-08-21 9/450 2026-08-22 11:28 by h4CP7TrQR8Lg
[基金申请] 2026国自然函评费到账 +9 羊腰板 2026-08-21 9/450 2026-08-22 00:03 by anji999
[基金申请] 放榜前的不淡定 40+4 snowwithsea 2026-08-19 14/700 2026-08-21 23:51 by cratir
[基金申请] 科研孤儿太难了 +17 我4大白菜 2026-08-20 18/900 2026-08-21 20:57 by zhangev
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 让我中一个面上吧! +11 大萍1987 2026-08-20 13/650 2026-08-21 14:48 by 20120902066
[基金申请] 感觉是下周放榜了 +7 angus9576 2026-08-17 12/600 2026-08-21 13:38 by weiyin
[基金申请] 我面上完蛋了 +7 且听虎啸 2026-08-20 8/400 2026-08-21 12:31 by 酷酷墨镜
[基金申请] 什么时候开奖? +7 CrisMessi 2026-08-18 7/350 2026-08-21 12:15 by weiyin
[基金申请] 今天基金会出结果吗?20260819 +15 kkkl_v 2026-08-19 16/800 2026-08-21 11:22 by 365372687
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[论文投稿] 投稿咨询 +5 wwm09 2026-08-17 7/350 2026-08-21 10:11 by 期刊论文帮手
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
[基金申请] 时间戳变了,能看出什么问题? +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[基金申请] 重要消息,中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
[基金申请] 今天放榜吗? +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
信息提示
请填处理意见