| 查看: 991 | 回复: 2 | ||||
apollosun9283金虫 (小有名气)
剑桥公爵
|
[求助]
常微分方程组初值问题的参数确定
|
|
方程如附件的图1和图2所示。 --------------------定义方程组--------------------------------------------------- function dydt=tangbin(t,y) global k1 k2 k3 k4 k5 k6 k7 f1=3*k1*y(1)-k2*y(1)*y(2)-k3*y(1)*y(3)-k4*y(1)*y(4)-... k5*y(1)*y(5)-k6*y(1)*y(6)-k7*y(1); f2=-1/3*k2*76/90*y(1)*y(2); f3=1/3*k2*106/90*y(1)*y(2)-1/3*k3*106/90*y(1)*y(3); f4=1/3*k3*136/90*y(1)*y(3)-1/3*k4*136/90*y(1)*y(4); f5=1/3*k4*166/90*y(1)*y(4)-1/3*k5*166/90*y(1)*y(5); f6=1/3*k5*196/90*y(1)*y(5)-1/3*k6*196/90*y(1)*y(6); f7=1/3*k6*226/90*y(1)*y(6)-1/3*k7*226/90*y(1); f8=1/3*k7*256/90*y(1); dydt=[f1;f2;f3;f4;f5;f6;f7;f8]; --------------------------以k值的猜测值计算y----------------------------- clear all;clc global k1 k2 k3 k4 k5 k6 k7 k1=5; k2=900; k3=500; k4=500; k5=500; k6=700; k7=3; y0=[0.25 0.75 0 0 0 0 0 0]; [t,y]=ode45(@tangbin,[1/30 1/24 1/18 1/12],y0); yexp=0.01*[2.21 62.89 19.04 6.86 2.34 0.80 0.27 0.02;... 1.82 59.76 20.23 8.07 2.97 1.02 0.33 0.05;... 1.15 55.64 22.27 9.48 3.73 1.34 0.47 0.12; 1.06 55.13 21.85 9.69 3.99 1.59 0.62 0.22]; s=sum(sum((y-yexp).^2)) ------------------------------------------------------------------------------- 现在我遇到的问题是，下一步如何调整k值，使得s=sum(sum((y-yexp).^2))尽可能的小。我设想过给k构建一个循环语句来调整k，但是如何调整，没有思路。 请高手帮帮忙！  1.jpg  2.jpg |
回复此楼» 猜你喜欢
号外号外
已经有4人回复
-
能否退出参与的面上项目解除限项
已经有22人回复
-
人气不行了
已经有9人回复
-
没有任何消息-是不是就凉了
已经有6人回复
-
如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》
已经有5人回复
-
2026国自然函评费到账
已经有19人回复
-
建议基金发布提前给出明确的时间点
已经有16人回复
-
今日不放榜？网传国自然预计 8 月 27 日可查结果
已经有21人回复
-
估计是周四
已经有4人回复
-
我面上完蛋了
已经有9人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab微分方程组参数拟合，以周为单位求解，汇总后以年为单位进行数值比较
已经有4人回复
- matlab数值求解边界条件微分方程组
已经有7人回复
- 高金求助matlab解微分方程组
已经有12人回复
- matlab解微分方程组
已经有15人回复
- LABVIEW四阶龙格库塔法求解一阶微分方程组
已经有12人回复
- 【求助】向各位大侠求助matlab求解微分方程组遇到的一个问题
已经有21人回复
- 【求助】请教非齐次常微分方程组的解析解法
已经有4人回复
- 【求助】求助一个一阶常微分方程的初值问题
已经有16人回复
- 【求助】急请微分方程高手看下这个方程组！
已经有14人回复
- 【求助】偏微分方程和常微分方程组成的方程组怎么求，大家帮帮忙啊
已经有3人回复
- 【求助】常微分方程组求解中系数与某变量值关联的问题
已经有12人回复
- 【求助】常微分方程组的解析解
已经有3人回复
- 【求助】如何用Runge-Kutta迭代求解二阶常微分方程组【已解决】
已经有9人回复
2楼2012-10-12 23:10:08
月只蓝
主管区长 (职业作家)
-
专家经验: +1059 - 应助: 1712 (讲师)
- 贵宾: 8.888
- 金币: 68171.7
- 散金: 1938
- 红花: 443
- 沙发: 4
- 帖子: 4373
- 在线: 3292.4小时
- 虫号: 1122189
- 注册: 2010-10-14
- 专业: 宇宙学
- 管辖: 计算模拟区
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
感谢参与，应助指数 +1
apollosun9283: 金币+150, ★★★★★最佳答案, 十分感谢！ 2012-10-13 09:06:40
xiegangmai: 金币+3 2012-10-14 08:23:29
感谢参与，应助指数 +1
apollosun9283: 金币+150, ★★★★★最佳答案, 十分感谢！ 2012-10-13 09:06:40
xiegangmai: 金币+3 2012-10-14 08:23:29
|
参考实用计算机模拟书，程序如下，做出来问题，至于效果的话，请楼主调一下k初值： function canshinihe clear all;clc tspan=[0 1/30 1/24 1/18 1/12]; y0=[0.25;0.75;0;0;0;0;0;0]; k0=[1 1 1 1 1 1 1]*0.1; lb=[0 0 0 0 0 0 0]; ub=[1 1 1 1 1 1 1]*10^5; yexp=0.01*[2.21 62.89 19.04 6.86 2.34 0.80 0.27 0.02;... 1.82 59.76 20.23 8.07 2.97 1.02 0.33 0.05;... 1.15 55.64 22.27 9.48 3.73 1.34 0.47 0.12; 1.06 55.13 21.85 9.69 3.99 1.59 0.62 0.22]; % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数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('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4)) fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5)) fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6)) fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7)) fprintf(' The sum of the squares is: %.1e\n\n',resnorm) %--------------------------------------------------------- function f = ObjFunc(k,tspan,y0,yexp) % 目标函数 [t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k); ysim(:,1) = Xsim(2:end,1); ysim(:,2) = Xsim(2:end,2); ysim(:,3) = Xsim(2:end,3); ysim(:,4) = Xsim(2:end,4); ysim(:,5) = Xsim(2:end,5); ysim(:,6) = Xsim(2:end,6); ysim(:,7) = Xsim(2:end,7); ysim(:,8) = Xsim(2:end,8); f =[ysim(:,1)-yexp(:,1); ysim(:,2)-yexp(:,2); ysim(:,3)-yexp(:,3);... ysim(:,4)-yexp(:,4); ysim(:,5)-yexp(:,5); ysim(:,6)-yexp(:,6);... ysim(:,7)-yexp(:,7); ysim(:,8)-yexp(:,8)]; %---------------------------------------------------------- function dydt = KineticsEqs(t,y,k) f1=3*k(1)*y(1)-k(2)*y(1)*y(2)-k(3)*y(1)*y(3)-k(4)*y(1)*y(4)-... k(5)*y(1)*y(5)-k(6)*y(1)*y(6)-k(7)*y(1); f2=-1/3*k(2)*76/90*y(1)*y(2); f3=1/3*k(2)*106/90*y(1)*y(2)-1/3*k(3)*106/90*y(1)*y(3); f4=1/3*k(3)*136/90*y(1)*y(3)-1/3*k(4)*136/90*y(1)*y(4); f5=1/3*k(4)*166/90*y(1)*y(4)-1/3*k(5)*166/90*y(1)*y(5); f6=1/3*k(5)*196/90*y(1)*y(5)-1/3*k(6)*196/90*y(1)*y(6); f7=1/3*k(6)*226/90*y(1)*y(6)-1/3*k(7)*226/90*y(1); f8=1/3*k(7)*256/90*y(1); dydt=[f1;f2;f3;f4;f5;f6;f7;f8]; |
3楼2012-10-13 07:13:09
