<DIV class=quote><B>以下是引用<I>irewine</I>在2006-4-28 10:22:03的发言:</B><br>感谢教授回答。<br>ODE45无法实现定步长对吧?<br>我现在毕业设计的一个关键步骤就是要龙格库塔法定步长解微分方程。有什么办法能在本机的matlab上使用ode4函数?</DIV><br>ode45可以实现等步长输出,但是计算的时候仍然是采用的是变步长<br><br><br>ode4.m<br><br>function [time, y_sol, y_drv] = ode4( dydt, time, y0, u, params )<br>% [time, y_sol, y_drv] = ode4( dydt, time, y0, u, params )<br>%<br>% Solve a system of nonhomogeneous ordinary differential equations using the <br>% 4th order Runge-Kutta method. <br>%<br>% Input Variable Description <br>% -------------- ----------- <br>% dydt : a function of the form y_dot = dydt(t,y,u,params)<br>% which provides the state derivative given <br>% the state, y, and the time, t <br>% time : a vector of points in time at which the solution <br>% is computed <br>% y0 : the initial value of the states <br>% params : vector of other parameters for dydt<br>% u : m-by-p matrix of forcing for the ode, where m<br>% m = number of inputs and n = length(time);<br>%<br>% Output Variable Description<br>% -------------- -----------<br>% time : is returned, un-changed<br>% y_sol : the solution to the differential equation at <br>% times given in the vector time.<br>%<br>% Henri Gavin, Civil and Environmental Engineering, Duke Univsersity, Jan. 2001<br>% WH Press, et al, Numerical Recipes in C, 1992, Section 16.1<br>points = length(time);<br>n = length(y0(:));<br>if nargin < 5, params = 0; end<br>if nargin < 4, u = zeros(1,points); end<br>y_sol = zeros(n,points); % memory allocation<br>y_drv = zeros(n,points); % memory allocation<br>y_sol(:,1) = y0(:); % the initial conditions<br>for p = 2:points<br>t = time(p); <br>dt = t - time(p-1); % the time step for this interval<br>dt2 = dt/2; % half of the time step<br>k1 = feval( dydt, t - dt, y0, u(:,p-1), params );<br>k2 = feval( dydt, t - dt2, y0 + k1*dt2, (u(:,p-1)+u(:,p))/2.0, params );<br>k3 = feval( dydt, t - dt2, y0 + k2*dt2, (u(:,p-1)+u(:,p))/2.0, params );<br>k4 = feval( dydt, t, y0 + k3*dt, u(:,p), params );<br>y0 = y0 + ( k1 + 2*(k2 + k3) + k4 )*dt/6; <br>y_sol(:,p) = y0(:);<br>y_drv(:,p-1) = k1(:);<br>end<br>
[此贴子已经被作者于2006-4-30 18:58:30编辑过]
|