Matlab讨论区 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

zhubo0709 发表于 2006-6-6 17:00

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

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

liuyh 发表于 2006-6-6 17:01

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

zhubo0709 发表于 2006-6-6 17:07

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

wangzi629 发表于 2006-6-6 17:17

是一阶常微分方程组吗？

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

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

zhubo0709 发表于 2006-6-6 19:24

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

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

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

happy 发表于 2006-6-6 21:58

回复：（zhubo0709）[求助]求个用4阶runge-kutta法解 ...

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

zhubo0709 发表于 2006-6-6 22:36

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

wangzi629 发表于 2006-6-7 09:39

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

grta 发表于 2006-6-7 10:06

是变步长的吗

zhubo0709 发表于 2006-6-7 13:41

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

wachen 发表于 2006-6-7 15:45

回复：（liuyh）我有用C++写的程序，需要的话，晚上可...

可不可以将你的c++程序传我一份，我在编c的杰微分方程组的程序，遇到了些麻烦，非常谢谢，qq:344244055..
[此贴子已经被作者于2006-6-8 9:00:15编辑过]

zhubo0709 发表于 2006-6-7 22:01

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

timchow 发表于 2007-1-11 16:16

没搞懂

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

packie 发表于 2007-1-14 10:42

改变一下tspan=;变成

页: [1]

声振论坛's Archiver

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

回复：（zhubo0709）[求助]求个用4阶runge-kutta法解 ...

回复：（liuyh）我有用C++写的程序，需要的话，晚上可...

没搞懂

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