如何用 matlab ode23s 求解 ODE 四阶？

Question

实际上我正在尝试使用 ode23s 求解以下 ODE。不幸的是我真的不知道该怎么做。

颂歌如下：

bo(t)*y + b1(t)*y' + b2(t)*y'' + b3(t)*y''' + b4(t)*y''' = ao(t)*cst

(那些导数应该用点标记，因为它们都是时间导数，例如 y' = dy/dt, y'' = d2y/dt2 等。参数 bo, b1, ..., b4 和 ao 作为 matlab 函数给出。'cst' 是常数。

首先，我将这个方程式转化为一个 ODE 系统，这导致我得到以下函数：

function [dy] = odetest(y,cst,ao,bo,b1,b2,b3,b4)
  dy = [ y(2); ...
         y(3); ...
         y(4); ...
         (-b1*y(2) + ao*cst - b2*y(3) - b3*y(4) + bo*y(1))*(1/(b4)) ];
end

所以...为了简化程序（并作为第一个近似值），我希望每个时间步长的时间相关参数 ao、bo、b1、...b4 保持不变。

现在问题开始了：我不知道如何将这个函数放在 ode23s 中，它可以在给定的初始条件下针对任意时间步长求解。也许我在时间导数的转换中犯了一个错误......

我猜想解决方案必须是这样的，不幸的是它不起作用

% file where the Ode has to be solved
% Loop over all time-steps

for i=1:number_timesteps

   time0 = i-1;
   time  = i;

   ao = function_ao(i); % function of parameter a, which is now constant in each timestep

   bo = function_bo(i);
   b1 = function_b1(i);
   b2 = function_b2(i);
   b3 = function_b3(i);
   b4 = function_b4(i);

   y = value1;   % initial conditions for y
   dy = value2;
   ddy = value3;
   dddy = value4;

   cst = value5; % value of cst in the actual timestep

   [T, Y] = ode23s(@odetest, [time0 time1] [y; dy; ddy; dddy], ...
   cst, ao, bo, b1, b2, b3, b4 ); % the arguments are the inital conditions, parameters ao, bo, ...

end

但是……我真的不明白我做错了什么。而且我也不知道如果参数 ao, bo, b1, ..., b4 在一个时间步内发生变化，所以 ode23s 必须以某种方式使用它们的功能。

如果有人能帮助我就太好了。提前致谢！

编辑 - 实施

我现在修改了我的代码如下，但是不幸的是报错了"Error using odetest (line 5) Not enough input arguments."不知道为什么:/

odetest.m:

    function [ dy ] = odetest( ttime, y, cst, ff1, ff2, ff3, ff4, ...
                  df1, df2, df3, df4, dd2, dd31, dd32, dd33, ...
                  dd34, ttk, tt0, EE1)

    EE0 = fe0(ff1, ff2, ff3, ff4, ttk, tt0, ttime);
    e1 = fe1(df1, dd2, dd31, ttk, tt0, ttime);
    e2 = fe2(df2, dd2, dd32, ttk, tt0, ttime);
    e3 = fe3(df3, dd2, dd33, ttk, tt0, ttime);
    e4 = fe4(df4, dd2, dd34, ttk, tt0, ttime);

    a0 = fa0(EE0, e1, e2, e3, e4);    
    b0 = fb0(e1, e2, e3, e4);
    b1 = fb1(EE0, EE1, EE1, EE1, EE1, e1, e2, e3, e4);
    b2 = fb2(EE0, EE1, EE1, EE1, EE1, e1, e2, e3, e4);
    b3 = fb3(EE0, EE1, EE1, EE1, EE1, e1, e2, e3, e4);
    b4 = fb4(EE0, EE1, EE1, EE1, EE1);

    dy = [ y(2); ...
           y(3); ...
           y(4); ...
           (-bb_1*y(2) + aa_0*cst - bb_2*y(3) - bb_3*y(4) + bb_0*y(1))*(1/(bb_4))];

end

文件本身，我调用这个函数的地方是

% something before (set constant parameters val_i and define initial
% conditions yy_t0 for first step ...
for i=1:n

   dbl_t0 = t(i-1,1);
   dbl_t  = t(i,1);

   cst = act(i,1) - act(i-1);


   [T, Y] = ode23s(@odetest, [dbl_t0 dbl_t], [yy_t0 ...
            dyy_t0 ddy_t0 dddy_t0],cst, ...
            val_f1, val_f2, val_f3, val_f4, val_d1, val_d2, val_d3, val_d4, val_dd2, ...
            val_dk1, val_dk2, val_dk3, val_dk4, val_tk, val_t0, val_E1);

   yy_t = Y(1);
   dyy_t = Y(2);
   ddyy_t = Y(3);
   dddyy_t = Y(4);

   yy_t0 = yy_t;
   dyy_t0 = dyy_t;
   ddyy_t0 = ddyy_t;
   dddyy_t0 = dddyy_t;

end

...有人看到这里有错误吗？

上次编辑

我发现了错误。 - 我必须按以下方式调用 ode23s：

[T, Y] = ode23s(@(ttime,y)odetest(ttime, y, cst, ...
            f_1, f_2, f_3, f_4, d_11, d_12, d_13, d_14, d_2, d_31, ...
            d_32, d_33, d_34, t_k, t_0, E_1), [dbl_t0 dbl_t], ...
            [yy_t0 dyy_t0 ddyy_t0 dddyy_t0]);

...现在可以使用了。

Answer 1

不要对单个时间步使用 ode23s，而是让它对整个时间序列进行积分，并利用您的 odetest 函数（微分系统函数）为每个时间步重新评估您的参数。您将获得更高的准确性，尽管它需要您自己进行工作以确保正确设置参数函数以适应您使用的时间跨度。对了，ODEFUN的第一个参数是时间，不是你的状态！

function [dy] = odetest(t, y, cst)
    ao = function_ao(t);
    b1 = function_b1(t);

    % ...re-evaluate all your parameters here...

    dy = [ y(2); ...
           y(3); ...
           y(4); ...
           (-b1*y(2) + ao*cst - b2*y(3) - b3*y(4) + bo*y(1))*(1/(b4)) ];
end

...

% Somewhere else in your code
y = value1;   % initial conditions for y
dy = value2;
ddy = value3;
dddy = value4;
[T, Y] = ode23s( @odetest, [0 number_of_timesteps], [y dy ddy dddy], cst );

这种风格在切换求解器时特别有用，因为它们会使用不同的积分技术来跳跃不同的时间位置（它们并不总是在特定时间点积分）。

如果您希望在 特定的 时间点进行积分，您还可以传递一个时间序列 ( 1:number_of_timesteps ) 而不是上面所示的时间跨度。

如何用 matlab ode23s 求解 ODE 四阶？

How to solve an ODE 4th order with matlab ode23s?

matlab

ode