matlab 一阶常微分方程求参数

方程是 (dx/dt)/x=r-(r/xm)*x
其中r和xm是两个常数参数,有一个表格数据已确定t跟x的关系(共17组数据)
要求使用这些数据确定r和xm.
详细说明matlab程序及大概原理。用到函数 ode45
另外要误差比较的函数
注:就一天时间,满意另加100分

正着做:给一个r和xm,通过ode45求得t,x,然后再与你的数据对比。最终选择一个合适的r和xm。就要一直变参数,我觉得比较难。至少来说r和xm的选择范围太大了。而且ode45得到的t,x与你的实验t,x肯定不是同一个t下面的数据,也不好比较。

不过基于你的方程比较简单,我们可以直接解出它的通解来:
>> xx=dsolve('Dx/x=r-r/b*x')
得到:
xx =
b/(1+exp(-r*t)*C1*b)

自己想手算的话,见附录。

===================================
===================================
以下拟合
然后再用你的17组数据,拟合出下面的r,b,c1来。
最终就可以确定你的r,xm了。这里的xm=b。
下面是通过数据来拟合
在Matlab下输入:edit,然后将下面两行百分号之间的内容,复制进去,保存
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=zhidao_joans321(a,t)
r=a(1);
b=a(2);
C1=a(3);

F=b./(1+exp(-r*t)*C1*b);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

你设好t,x后,比如
t=1:17;
x=sin(t);
%在Matlab下输入:

[A,res]=lsqcurvefit('zhidao_joans321',ones(1,3),t,x);
A
r=A(1)
xm=A(2)

================================
================================
附录:
xm我用b代替
(dx/dt)/x/(r-r/b*x)=1
dx/[rx(1-1/b*x)]=dt
[1/x+(1/b)/(1-x/b)]dx=r*dt
积分得
ln(x)-ln(1-x/b)=r*t+C
-ln[(1-x/b)/x]=r*t+C
ln[1/x-1/b]=-r*t-C
两边指数得
1/x-1/b=C1*exp(-r*t)
1/x=1/b+C1*exp(-r*t)
x=b/[1+b*C1*exp(-r*t)]

x=b/(1+exp(-r*t)*C1*b)
温馨提示:答案为网友推荐,仅供参考
第1个回答  2015-11-07

 ‘的含义为转置,*为乘法;n'*k,就代表矩阵n转置后与矩阵k相乘


参考一:计算离散时间傅里叶变换,并绘制图形。
已知有限长序列x(n)={1,2,3,4,5}。
n=-1:3;

x=1:5;

k=0:500;

w=(pi/500)*k;

X=x*(exp(-j*2*pi/500)).^(n'*k);

subplot(211)

plot(w/pi,imag(X))

 subplot(212)

plot(w/pi,abs(X))

grid on

向左转|向右转
第2个回答  2015-10-07
给一个r和xm,通过ode45求得t,x,然后再与你的数据对比。最终选择一个合适的r和xm。就要一直变参数,我觉得比较难。至少来说r和xm的选择范围太大了。而且ode45得到的t,x与你的实验t,x肯定不是同一个t下面的数据,也不好比较。
不过基于你的方程比较简单,我们可以直接解出它的通解来:
>> xx=dsolve('Dx/x=r-r/b*x')
得到:
xx =
b/(1+exp(-r*t)*C1*b)
自己想手算的话,见附录。
===================================
===================================
以下拟合
然后再用你的17组数据,拟合出下面的r,b,c1来。
最终就可以确定你的r,xm了。这里的xm=b。
下面是通过数据来拟合
在Matlab下输入:edit,然后将下面两行百分号之间的内容,复制进去,保存
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=zhidao_joans321(a,t)
r=a(1);
b=a(2);
C1=a(3);
F=b./(1+exp(-r*t)*C1*b);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
你设好t,x后,比如
t=1:17;
x=sin(t);
%在Matlab下输入:
[A,res]=lsqcurvefit('zhidao_joans321',ones(1,3),t,x);
A
r=A(1)
xm=A(2)
================================
================================
第3个回答  2008-09-24
他做的是什么~~~