&& 查看话题
用matlab解一个含参数的一元三次方程，高手帮忙看看这个程序有什么问题，谢谢！
用matlab解一个含参数的一元三次方程：
C3*P(z)^3+C1*P(z)+D(z)-P_av=0&&
z是一线性数组，C3，C1常数，P_av是P在不同z的根的平均值。
思路是首先假设一个初始值P_av0, 然后解得P的一组根，取平均值和P_av0比较，如果收敛，则为所求的根，否则把新的平均值代回方程，再求根，直到收敛。
程序如下：
function P_real=P_distribution(a,b,e,q,d,N,m,z,P_r,Dsc)
format long
% define parameters
a=-6.36e7;& && && && && && && && && && && && && && && && && &
b=1.38e9;& && && && && && && && && && && && && && && && && &
e=8.85e-11;& && && && && && && && && && && && && && && && && && &
q=1.61e-19;& && && && && && && && && && && && && && && && &&&
d=1e-7;& && && && && && && && && && && && && && && && && &&&
% coeffecients
% start value
P_av0=0.2;& && && & % initial value for the P average
P_real=zeros(1,11);
%start the caculation
& &while 1
& && & z=m*1e-8;
& && & Dsc=-q*N*(d/2-z);
& && & c0=Dsc-P_av0;
& && & P_eq=;
& && & P_r=roots(P_eq);& &%get the roots
& && & for k=1:length(P_r)
& && && &&&P_rr=isreal(P_r(k));& &%get the real roots
& && & end
& && & pp=P_r(P_rr);
& && & P_real(1,m+1)=abs(pp);& && & %assign the real roots to the matrix
& && & m=m+1;& && && &&&%solve the equation again for another value
& && & if m>10
& && && &&&break
& && & end
& && & P_av1=mean(P_real);& & %get the average value
& && & if abs(P_av1-P_av0)<1e-11& &%less than the tolerance
& && && &&&break
& && & else
& && && &&&P_av0=P_av1;
& && & end
程序在matlab中无法运行，不知道问题出在哪里，高手帮忙看看，多谢了！
if abs(P_av1-P_av0)<1e-11& &%less than the tolerance
中精度1e-11设置的太高了，关键与前面的roots
P_r=roots(P_eq);& &%get the roots
相比精度太高了，
建议把root换成fzero，因为fzero可以直接控制误差的 if m>10
end : Originally posted by wangww2011 at
if abs(P_av1-P_av0)&1e-11& &%less than the tolerance
中精度1e-11设置的太高了，关键与前面的roots
P_r=roots(P_eq);& &%get the roots
相比精度太高了，
建议把root换成fzero，因为fzero可以直接 ... 尝试了降低精度，和用fzero，
但是还是无法在matlab中运行，matlab也没有报错，不知道问题出在哪里 : Originally posted by csgt0 at
end m>10已经给出条件了，
所以并不需要m=0呀 : Originally posted by xmuxiaoyu at
m&10已经给出条件了，
所以并不需要m=0呀... 1.是可以运行的
2.是因为你原意是将每次preal赋给下一次的P_av1，但是你的程序是把所有的preal均值赋给pav，因为你的m一直增大，preal就越来越大。降低精度就可以看到很快就算完了。 : Originally posted by xmuxiaoyu at
m&10已经给出条件了，
所以并不需要m=0呀... 楼上说的对&&你的m是一直增大的 应该加上m=0才行的
你可以在两个while循环之间加上m=0，不需要改精度 : Originally posted by csgt0 at
1.是可以运行的
2.是因为你原意是将每次preal赋给下一次的P_av1，但是你的程序是把所有的preal均值赋给pav，因为你的m一直增大，preal就越来越大。降低精度就可以看到很快就算完了。... 我的意思是把P_real数组的平均值赋给P_av1。
m增大，只是求不同z值的根，也就是一组方程的根。一共11个值。
P_real不会越来越大的，只是给数组里的每一个数赋予一个根。
或者我没有理解你的意思，这个程序能在你的机上运行吗？十分感谢！ : Originally posted by xmuxiaoyu at
我的意思是把P_real数组的平均值赋给P_av1。
m增大，只是求不同z值的根，也就是一组方程的根。一共11个值。
P_real不会越来越大的，只是给数组里的每一个数赋予一个根。
或者我没有理解你的意思，这个程序能在 ... 精度降到E-51改不改m都能运行，但你不加m=0跟你的愿意是不同的 降低精度都能运行 : Originally posted by xmuxiaoyu at
我的意思是把P_real数组的平均值赋给P_av1。
m增大，只是求不同z值的根，也就是一组方程的根。一共11个值。
P_real不会越来越大的，只是给数组里的每一个数赋予一个根。
或者我没有理解你的意思，这个程序能在 ... 你的程序不加上m=0根本就不是你想要的，很明确的说就是错的，虽然降低精度会出现结果，但结果也是错的，所以还是请楼主加上m=0。 : Originally posted by csgt0 at
精度降到E-51改不改m都能运行，但你不加m=0跟你的愿意是不同的... 加了m=0 可以计出结果，十分感谢！
但是这个结果跟之前在 Igor Pro 中计得的不太一样，
所以还是有点疑惑，这个程序是不是计算我要的东西。 : Originally posted by wangww2011 at
你的程序不加上m=0根本就不是你想要的，很明确的说就是错的，虽然降低精度会出现结果，但结果也是错的，所以还是请楼主加上m=0。... 多谢提醒！
现在能够得到结果。
还在考虑这个结果是不是我要算的东西。
var cpro_id = 'u1216994';
欢迎监督和反馈：本帖内容由
提供，小木虫仅提供交流平台，不对该内容负责。欢迎协助我们监督管理，共同维护互联网健康，如果您对该内容有异议，请立即发邮件到
联系通知管理员，也可以通过QQ周知，我们的QQ号为：8835100
我们保证在1个工作日内给予处理和答复，谢谢您的监督。
小木虫，学术科研第一站，为中国学术科研研究提供免费动力
广告投放请联系QQ： &
违规贴举报删除请联系邮箱： 或者 QQ:8835100
Copyright &
eMuch.net, All Rights Reserved. 小木虫 版权所有&& 查看话题
一种高质量拟合常微分方程参数的方法：低版本1stOpt和MATLAB联用
就方程参数拟合问题而言，常见非线性方程（组）参数拟合和常微分方程（组）参数拟合两种。
本帖主要讨论常微分方程参数的拟合问题。限于本人水平，存在纰漏甚至错误在所难免，欢迎大家批评指正，同时希望能抛砖引玉。
1.影响拟合质量的因素
Origin软件的Help文档总结了影响拟合质量的主要因素，总共5条，列出如下：
（1）参数初值
（2）模型（公式）的正确性或者适用性
（3）迭代计算的中止的判断指标高低
（4）数据点的充分性
（5）数据点的噪声
实际上，拟合算法的先进性同样影响拟合的质量，在此基础上，加上一条：
（6）算法先进与否及其适用性
本帖主要讨论第一个因素，即参数初值，对拟合质量的影响。在本体后续的讨论中可以看到，参数初值的选取，对最终拟合结果影响相当显著。
2.本帖讨论的例子出处
本帖讨论的例子出处见：http://emuch.net/bbs/viewthread.php?tid=7595881&fpage=1
直观起见，对上述例子中的问题表述如下：
已知自变量t，因变量y，y和t满足一下函数式：
dy/dt=0./(4.*k1*(1-0.02119*y)^k2)+1/k3/45727)；（需要指出的是，该公式中参数可进一步合并，在本帖中并未作合并处理）
拟合出参数k1，k2，k3。
3.结果与讨论
运行MATLAB软件，调用lsqnonlin函数和ode45函数，参数初值按k1=k2=k3=1选取，程序代码如下：
function KineticsEst1_Int_modified_by_Yuezhilan
format long
tspan=-31;
k0=;& &%%%请注意这里，初值的选取
y0=20.4122;
& & lsqnonlin(@ObjFunc,k0,lb,ub,,tspan,y0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t待拟合参数 k1 = %.6f\n',k(1))
fprintf('\t待拟合参数 k2 = %.6f\n',k(2))
fprintf('\t待拟合参数 k3 = %.6f\n',k(3))
fprintf(' \t残差平方和= %.6f\n\n',resnorm)
ts=0:1:max(tspan);
=ode45(@KineticsEqs,ts,y0,,k);
= ode45(@KineticsEqs,tspan,y0,,k);
y=XXsim(2:end);
R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);
fprintf('\n\t决定系数R-Square = %.6f',R2);
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best');
plot(tspan(2:end),yr,'r*',,),axis();
plot(yexp,y,'ro',,,'b-');
%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)& && && &&&
= ode45(@KineticsEqs,tspan,y0,,k) ;
ysim = Xsim(2:end);
size(ysim);
size(yexp);
%----------------------------------------------------------
function dydt = KineticsEqs(t,y,k)
beta(1)=k(1);
beta(2)=k(2);
beta(3)=k(3);
dydt = 0./(4.*beta(1)*(1-0.02119*y)^beta(2))+1/beta(3)/45727);
计算结果：
使用函数lsqnonlin()估计得到的参数值为:
& & & & 待拟合参数 k1 = 585.320105
& & & & 待拟合参数 k2 = 14.003057
& & & & 待拟合参数 k3 = 0.065162
& & & & 残差平方和= 1.705667
& & & & 决定系数R-Square = 0.945774
Warning: Imaginary parts of complex X and/or Y arguments ignored
& In KineticsEst1_Int_modified_by_Yuezhilan at 30
Warning: Imaginary parts of complex X and/or Y arguments ignored
& In KineticsEst1_Int_modified_by_Yuezhilan at 33
Warning: Imaginary parts of complex X and/or Y arguments ignored
& In KineticsEst1_Int_modified_by_Yuezhilan at 35
由结果可见，决定系数R-Square = 0.945774，拟合质量一般，同时MATLAB警告结果含有复数。
通常而言，在作拟合的时候，常赋予初值诸如（0 0 0）、（1 1 1）、（0 1 1）……之类比较简单的数值以试探；更合理的方法，是事先对数据和公式进行观察，猜测初值的取值范围。不过，对于拟合公式较复杂的情况，通过观察实验数据很难判断初值的范围，比如在本例中，通过观察y随t的变化趋势可见（见附图1），基本上y随t递增，也就是说，dy/dt在t取值范围内应&0，仅能大概判断k1、k2和k3大于0的可能性比较大。
继续尝试其它初值的选取，结果列出如下
（k1 k2 k3）& && && && && && && && && && && && && && && & 决定系数R-square
（1 0 0）& && && && && && && && && && && && && && && && && &0.958200
（0 1 0）& && && && && && && && && && && && && && && && && &0.000000（拟合失败）
（0 0 1）& && && && && && && && && && && && && && && && && &0.957533
（1 1 0）& && && && && && && && && && && && && && && && && &0.943656
（1 0 1）& && && && && && && && && && && && && && && && && &0.962539
（0.5 1 1）& && && && && && && && && && && && && && && && &0.918617
（1 0.5 1）& && && && && && && && && && && && && && && && &0.958019
（1 1 0.5）& && && && && && && && && && && && && && && && &0.919118
（0.5 0.5 0.5）& && && && && && && && && && && && && && &0.955740
& &……& && && && && && && && && && && && && && && && && && && &&&……
由此可见，不同的初值取法，拟合结果差异显著，且在上述尝试的结果中，决定系数最高仅为0.962539。
1stOpt软件处理拟合问题强大高效，使用者在操作层面上不需要赋予初值，即可获得极佳的拟合效果。低版本的1stOpt（比如1.5）软件并不支持常微分方程的拟合，但其对非线性方程的拟合已经足够强大。
为了利用低版本1stOpt软件对代数方程强大的拟合能力，对本例问题初值的选取提出了一个解决方法，步骤如下：
1.MATLAB计算出数据的数值导数dy/dt；
2.记dy/dt=Y，原常微分方程化为Y=f（t，y，k1，k2，k3）的二元代数方程，1stOpt拟合出各个k值；
3.将上述拟合出的k值作为初值，输入到上述MATLAB程序中。
上述步骤1中求数值导数的MATLAB程序如下：
function initial_value_determination
sp=spap2(knots,K,t,yexp);
pp=fnder(sp);
dydt=fnval(pp,t);
计算结果如下：
& && && && && && & 0& &2.042&&20.999
& &1.000& &1.365&&22.999
& &2.000& &1.689&&24.999
& &3.000& &1.012&&25.001
& &4.000& &0.335&&26.001
& &5.000& &0.659&&27.002
& &6.000& &0.037&&27.998
& &7.000& &0.469&&27.001
& &8.000& &0.901&&27.998
& &9.000& &0.332&&27.000
&&10.000& &0.350&&27.999
&&11.000& &0.367&&27.000
&&12.000& &0.384&&27.000
&&13.000& &0.402&&27.000
&&14.000& &0.419&&27.998
上述步骤2中1stOpt代码和计算结果如下：
Parameters k1,k2,k3;
Variable& &t,dydt,y;
Function& &dydt = 0./(4.*k1*(1-0.02119*y)^k2)+1/k3/45727);
& && && && && && & 0& &2.042&&20.999
& &1.000& &1.365&&22.999
& &2.000& &1.689&&24.999
& &3.000& &1.012&&25.001
& &4.000& &0.335&&26.001
& &5.000& &0.659&&27.002
& &6.000& &0.037&&27.998
& &7.000& &0.469&&27.001
& &8.000& &0.901&&27.998
& &9.000& &0.332&&27.000
&&10.000& &0.350&&27.999
&&11.000& &0.367&&27.000
&&12.000& &0.384&&27.000
&&13.000& &0.402&&27.000
&&14.000& &0.419&&27.998
均方差(RMSE): 0.567
残差平方和(SSE): 0.613
相关系数(R): 0.64
相关系数之平方(R^2): 0.779
决定系数(DC): 0.033
卡方系数(Chi-Square): 0.371
F统计(F-Statistic): 517.
参数& & & & 最佳估算
----------& & & & -------------
k1& && && && & & & 9094
k2& && && && & & & 26.4
k3& && && && & & & 0.0055
将上述获得的k值输入，计算结果为：
使用函数lsqnonlin()估计得到的参数值为:
& & & & 待拟合参数 k1 = 7508
& & & & 待拟合参数 k2 = 30.541351
& & & & 待拟合参数 k3 = 0.033981
& & & & 残差平方和= 0.177590
& & & & 决定系数R-Square = 0.994276&&
拟合效果见附图2、附图3、附图4。
由结果可见，决定系数R-Square = 0.994276，相比于之前有显著提高，可达到两个“9”的拟合质量。
祝福，祈福，愿楼主今后一切顺利，愿自己及女朋友工作顺利，祈福！ 愿楼主心想事成
var cpro_id = 'u1216994';
欢迎监督和反馈：本帖内容由
提供，小木虫仅提供交流平台，不对该内容负责。欢迎协助我们监督管理，共同维护互联网健康，如果您对该内容有异议，请立即发邮件到
联系通知管理员，也可以通过QQ周知，我们的QQ号为：8835100
我们保证在1个工作日内给予处理和答复，谢谢您的监督。
小木虫，学术科研第一站，为中国学术科研研究提供免费动力
广告投放请联系QQ： &
违规贴举报删除请联系邮箱： 或者 QQ:8835100
Copyright &
eMuch.net, All Rights Reserved. 小木虫 版权所有各位MATLAB高手，有一用MATLAB求解非线性方程组的程序题不会，希望帮帮小妹，解燃眉之急！！！_百度知道
各位MATLAB高手，有一用MATLAB求解非线性方程组的程序题不会，希望帮帮小妹，解燃眉之急！！！
求出120组r：N(t)=K&#47，每两个已知数带入方程;(1+(K&#47：54 65 50 48 43 48 41 44 48 51 44 46 43 45 42 39，N(1)至N(16)这16个已知数。求以N0为基数，并求出r;N0-1)*exp(-r*t))（该式子在此用MATLAB语言给出）已知：N0=43：给出一个式子、k的值、k的平均数题目,N(1)-N(16)的值分别为
是从1到16，即N(1)-N(16)都给出来了，不过N(t)的值已知K是式子中的一个未知数，能求出120组K和r。t已知，N(t)也是式子中的一个数，然后再求r和K的平均值。要求的是两个未知参数K和r。这个题目一共有120组式子？可以吗。这样说清楚点儿了吗
43-1)*exp(-r*(t+m)))&#39。K&#39;43-1)*exp(-r*t))&#39:16
m=1,K]=solve(&#39,&#39;;N=[54 65 50 48 43 48 41 44 48 51 44 46 43 45 42 39];r&#39还是有点问题,&#39;N(t)=K/for t=1？我一时改不过来,&#39，你也改一下;;:16-t
[r;(1+(K/(1+(K/N(t+m)=K&#47，你看看这个程序行吗
其他类似问题
其他2条回答
用插值法做呢
我好好想想。。。。
等待您来回答
下载知道APP
随时随地咨询
出门在外也不愁我在学习MATLAB的使用，想用有限差分的方法求解二阶偏微分方程，请高手指点！
[问题点数：100分]
我在学习MATLAB的使用，想用有限差分的方法求解二阶偏微分方程，请高手指点！
[问题点数：100分]
不显示删除回复
显示所有回复
显示星级回复
显示得分回复
只显示楼主
相关帖子推荐：
本帖子已过去太久远了，不再提供回复功能。&& 查看话题
求助，用matlab解微分方程组，希望高手能给予指点
微分方程的形式如下所示：
d2y/dx2 =a*sinh(y)
边界条件为y(0)=b,dy/dx(x=0处)=c。
求y与x的曲线关系。
因为已经知道y 不存在解析解，只存在数值解，而且文献中提到是用matlab中的bvp4c路径做的，所以请matlab高手和数学高手驻足，帮忙解决一下，不胜感激。
用1stopt应该也是可以求解的~ /view/953dce.html
/view/690deff95436.html
可以参考一下上面两个连接 : Originally posted by xxz903 at
用1stopt应该也是可以求解的~ 具体怎么做，可不可以指教？ : Originally posted by arising2010 at
/view/953dce.html
/view/690deff95436.html
可以参考一下上面两个连接 谢谢，我以前就看过了，可是对应在我的方程中时，就不知道是什么了 你的方程少个边界条件吧，是不是少个&dy/dx(x=某数值)= ”的边界条件？ 不然按你目前的边界条件，都是x=0处的初始值，这个微分方程已经不是边值问题了，而是初值问题了。在现有条件下，解该问题用ode可解，程序如下，复制进一个m文件，F5运行即可，根据需要输入a，b，c的数值。a=b=c=1的结果见附图1，该方程组有刚性，不知道你希望x什么多大的区间内求解。补充完边界条件后，可再发上来。
function solvebvp333
global a b c
a=input('参数a=');
b=input('参数b=');
c=input('参数c=');
xspan=0:0.02:1;
=ode15s(@fun332,xspan,y0);
plot(x,y(:,1),'bo-');
function dydx=fun332(x,y)
附图1.jpg : Originally posted by musejianglin at
具体怎么做，可不可以指教？... 1stopt是一款软件，3.0以上版本的有求解常微分方程初边值问题的功能，这款软件是需要购买的。 不必用1stopt， MATLAB 完全可以做到
这里是程序 ：
===================================
% main program
t0& & = 0;& &tend=5000;
tspan = t0:0.5:
v_y_ini = ;
= ode45(@ode_sys, tspan, v_y_ini);
v = Y(:, 1);
y = Y(:, 2);
===================================
% ode system
function&&= ode_sys( ~, Y )
v = Y(1);& &y = Y(2);
dy = zeros(2,1);
dy(1)&&= a * sinh( y );
=================================== 记得把
function&&= ode_sys( ~, Y )
单独另存为一个M文件 : Originally posted by 月只蓝 at
你的方程少个边界条件吧，是不是少个&dy/dx(x=某数值)= ”的边界条件？ 不然按你目前的边界条件，都是x=0处的初始值，这个微分方程已经不是边值问题了，而是初值问题了。在现有条件下，解该问题用ode可解，程序 ... 谢谢，高手。 : Originally posted by 月只蓝 at
你的方程少个边界条件吧，是不是少个&dy/dx(x=某数值)= ”的边界条件？ 不然按你目前的边界条件，都是x=0处的初始值，这个微分方程已经不是边值问题了，而是初值问题了。在现有条件下，解该问题用ode可解，程序 ... 谢谢，高手。能再帮我看看么？能不能弄出来？
附件.png : Originally posted by somomo91 at
不必用1stopt， MATLAB 完全可以做到
这里是程序 ：
===================================
% main program
t0& & = 0;& &tend=5000;
tspan = t0:0.5:
... 谢谢你了，那能不能再帮我弄弄这个。谢谢了
附件.png : Originally posted by musejianglin at
谢谢，高手。能再帮我看看么？能不能弄出来？
... 边界条件还是有问题吧，边界条件顾名思义是将解确定在一定的区间之内，你的边界条件只有右端，没有左端，怎么确定解呢，这样子说不定有无穷多个解。
能不能给出比如&y(r=b)=& 或者&dy/dr(r=b)=&(b&a)呢？ : Originally posted by 月只蓝 at
边界条件还是有问题吧，边界条件顾名思义是将解确定在一定的区间之内，你的边界条件只有右端，没有左端，怎么确定解呢，这样子说不定有无穷多个解。
能不能给出比如&y(r=b)=& 或者&dy/dr(r=b)=&qu ... 哦，这是个球，前面的条件是在r=a处即球的表面处，而在0&r&a范围内都是0.这样可以了么？ : Originally posted by musejianglin at
哦，这是个球，前面的条件是在r=a处即球的表面处，而在0&r&a范围内都是0.这样可以了么？... 是y(r=0)=0，还是dy/dr(r=0)=0呢？
不可能说0&r&a范围内都是0，如果这样的话，这个区间内原函数的数值全部确定了，都不需要解方程了。 : Originally posted by musejianglin at
哦，这是个球，前面的条件是在r=a处即球的表面处，而在0&r&a范围内都是0.这样可以了么？... 你的原始模型的方程一定要建对，不然解出来的结果都是没有意义的。 : Originally posted by 月只蓝 at
是y(r=0)=0，还是dy/dr(r=0)=0呢？
不可能说0&r&a范围内都是0，如果这样的话，这个区间内原函数的数值全部确定了，都不需要解方程了。... 是 y(r=0)=0. : Originally posted by 月只蓝 at
你的原始模型的方程一定要建对，不然解出来的结果都是没有意义的。... 是 y(r=0)=0
var cpro_id = 'u1216994';
欢迎监督和反馈：本帖内容由
提供，小木虫仅提供交流平台，不对该内容负责。欢迎协助我们监督管理，共同维护互联网健康，如果您对该内容有异议，请立即发邮件到
联系通知管理员，也可以通过QQ周知，我们的QQ号为：8835100
我们保证在1个工作日内给予处理和答复，谢谢您的监督。
小木虫，学术科研第一站，为中国学术科研研究提供免费动力
广告投放请联系QQ： &
违规贴举报删除请联系邮箱： 或者 QQ:8835100
Copyright &
eMuch.net, All Rights Reserved. 小木虫 版权所有}