[求助]求个用4阶runge-kutta法解 微分方程组的程序。来大虾帮帮忙！

zhubo0709

<P>请大家帮我一下。需要能解3个变量，也就是3个方程的程序。谢谢啦！！！！<BR><BR><BR>PS：不是用MATLAB自带的ODE函数。: )</P>

liuyh

<P>我有用C++写的程序，需要的话，晚上可以回去给你找</P>

zhubo0709

<P>C++我还不太熟！可能看不懂哦！不过谢谢2楼的好意啦！</P>

wangzi629

是一阶常微分方程组吗？

我以前发过一个解一阶常微分方程组的帖子，你可以去看看，看能用上不。

[ 本帖最后由 ChaChing 于 2009-4-17 08:23 编辑 ]

zhubo0709

回4楼：是一阶常微分方程组！
帖子我去找找看！

你的帖子只是一个方程的程序,不是方程组的!

[ 本帖最后由 ChaChing 于 2009-4-17 08:26 编辑 ]

happy

<P>function Y = ode4(odefun,tspan,y0,varargin)<BR>%ODE4  Solve differential equations with a non-adaptive method of order 4.<BR>%   Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates <BR>%   the system of differential equations y' = f(t,y) by stepping from T0 to <BR>%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.<BR>%   The vector Y0 is the initial conditions at T0. Each row in the solution <BR>%   array Y corresponds to a time specified in TSPAN.<BR>%<BR>%   Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters <BR>%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). <BR>%<BR>%   This is a non-adaptive solver. The step sequence is determined by TSPAN<BR>%   but the derivative function ODEFUN is evaluated multiple times per step.<BR>%   The solver implements the classical Runge-Kutta method of order 4.   <BR>%<BR>%   Example <BR>%         tspan = 0:0.1:20;<BR>%         y = ode4(@vdp1,tspan,[2 0]);  <BR>%         plot(tspan,y(:,1));<BR>%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, <BR>%     and plots the first component of the solution.   <BR>%</P>
<P>if ~isnumeric(tspan)<BR>  error('TSPAN should be a vector of integration steps.');<BR>end</P>
<P>if ~isnumeric(y0)<BR>  error('Y0 should be a vector of initial conditions.');<BR>end</P>
<P>h = diff(tspan);<BR>if any(sign(h(1))*h &lt;= 0)<BR>  error('Entries of TSPAN are not in order.') <BR>end  </P>
<P>try<BR>  f0 = feval(odefun,tspan(1),y0,varargin{:});<BR>catch<BR>  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];<BR>  error(msg);  <BR>end  </P>
<P>y0 = y0(:);   % Make a column vector.<BR>if ~isequal(size(y0),size(f0))<BR>  error('Inconsistent sizes of Y0 and f(t0,y0).');<BR>end  </P>
<P>neq = length(y0);<BR>N = length(tspan);<BR>Y = zeros(neq,N);<BR>F = zeros(neq,4);</P>
<P>Y(:,1) = y0;<BR>for i = 2:N<BR>  ti = tspan(i-1); if ~mod(i,100), disp(['t = ' num2str(ti)]), end<BR>  hi = h(i-1);<BR>  yi = Y(:,i-1);<BR>  F(:,1) = feval(odefun,ti,yi,varargin{:});<BR>  F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});<BR>  F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});  <BR>  F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});<BR>  Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));<BR>end<BR>Y = Y.';<BR></P>

zhubo0709

教授果然厉害啊! <BR>我区研究一下!有不懂的地方,再来问问大家!<BR><BR>太感谢拉!: )

wangzi629

在我的方程中x，y，fx等可以是向量。比如都取为x，y都取为3维列向量，这样该向量的每个元素就是一个变量。就是3个方程。

grta

是变步长的吗

zhubo0709

function f=CY(t,y)<BR><BR>V=y(1);n=y(2);p=y(3);<BR><BR>am=0.1*(25 + V)/(1 - exp(-0.1*V - 2.5));<BR>bm=4*exp(-(V + 50)/18);<BR>ah=0.07*exp(-0.05*V - 2.5);<BR>bh=1/(1 + exp(-0.1*V - 2));<BR>an=0.01*(20 + V)/(1 - exp(-0.1*V - 2));<BR>bn=0.125*exp(-(V + 30)/80);<BR>gI=1800;   gK=1700;   gP=11.5;<BR>gL=7;      VI=100;    VK= -75;<BR>VL=-40;    KC=0.18;   tx_n=0.00435;<BR>tp=5.0;<BR>minf=am/(am + bm);hinf=ah/(ah + bh);<BR>ninf=an/(an + bn);t_n=tx_n/(an + bn);<BR>   <BR>dVdt=gI*(minf^3)*(hinf)*(VI- V) + gK*(n^4)*(VK- V) + gP*p*(VK- V) + gL*(VL- V);<BR>dndt=(ninf - n)/t_n;<BR>dpdt=(((minf^3)*(hinf)*(VI - V) - KC*p/(1-p))/tp)*((1-p).^2);<BR><BR>f = [dVdt dndt dpdt ]';<BR><BR>上面是我方程的函数了！不知道能用9楼的办法不能?方程确实有点复杂了，有点搞不定了！-.-!希望大家赐教.

wachen

可不可以将你的c++程序传我一份，我在编c的杰微分方程组的程序，遇到了些麻烦，非常谢谢，qq:344244055..

[此贴子已经被作者于2006-6-8 9:00:15编辑过]

zhubo0709

function   ODDD<br>y0=[-45,0.12,0.45];<br>tspan=[0:1:100];<br>N = length(tspan);<br>h = ones(1,N);<br>y0 = y0(:);<br>neq = length(y0);<br>Y = zeros(neq,N);<br>F = zeros(neq,4);<br>Y（:,1)=y0;<br>for i=2:N<br>ti=tspan(i-1);<br>hi=h(i-1);<br>yi=Y(:,i-1);<br>F(:,1)=feval(@odefun,ti,yi);<br>F(:,2)=feval(@odefun,ti+0.5*hi,yi+0.5*hi*F(:,1));<br>F(:,3)=feval(@odefun,ti+0.5*hi,yi+0.5*hi*F(:,2));<br>F(:,4)=feval(@odefun,ti,yi+hi*F(:,3));<br>Y(:,i)=yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));<br>end<br>Y = Y.'<br><br>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>function f=odefun(t,s)<br>V=s(1,1);<br>n=s(2,1);<br>p=s(3,1);<br>am=0.1*(25 + V)/(1 - exp(-0.1*V - 2.5));<br>bm=4*exp(-(V + 50)/18);<br>ah=0.07*exp(-0.05*V - 2.5);<br>bh=1/(1 + exp(-0.1*V - 2));<br>an=0.01*(20 + V)/(1 - exp(-0.1*V - 2));<br>bn=0.125*exp(-(V + 30)/80);<br>gI=1800;   gK=1700;   gP=11.5;<br>gL=7;      VI=100;    VK= -75;<br>VL=-40;    KC=0.18;   tx_n=0.00435;<br>tp=5.0;<br>minf=am/(am + bm);hinf=ah/(ah + bh);<br>ninf=an/(an + bn);t_n=tx_n/(an + bn);   <br>dVdt=gI*(minf^3)*(hinf)*(VI- V) + gK*(n^4)*(VK- V) + gP*p*(VK- V) + gL*(VL- V);<br>dndt=(ninf - n)/t_n;<br>dpdt=(((minf^3)*(hinf)*(VI - V) - KC*p/(1-p))/tp)*((1-p).^2);<br>f = [dVdt dndt dpdt ]';<br>这个是我根据给的模型遍的程序了！可是运行后有警告，而且结果出不来。大虾帮我看看程序。<br><br>运行后的警告是：<br>&gt;&gt; Warning: Divide by zero.<br>(Type "warning off MATLAB:divideByZero" to suppress this warning.)<br>&gt; In C:\MATLAB6p5\work\ODDD.m (odefun) at line 41<br>  In C:\MATLAB6p5\work\ODDD.m at line 17<br><br>

[此贴子已经被作者于2006-6-7 22:08:20编辑过]

timchow

楼主问题解决了没有
可否再讨论下

packie

改变一下tspan=[0:1:100];变成[0.1：1：100]