0

我在 t=0 时求解一组 ODE (dy/dt)，所有初始条件 t=0 y_0=(0,0,0)。我可以在不同时间向 y 值添加一些数字（例如，在 t=10 时，应将 y1 添加到该数字；在 t=20 时，应将 y2 添加到该数字等）并求解方程？

4

2 回答 2

2

以您建议的方式（以及@macduff 说明的方式）在您的 ODE 中插入较大的不连续性可能会导致精度降低和计算时间延长（尤其是ode45-ode15s可能是更好的选择，或者至少确保您的绝对和相对公差是合适的）。你有效地产生了一个非常僵硬的系统。如果您想从特定时间开始向 ODE 添加一些数字，请记住求解器仅在其自己选择的特定时间点评估这些方程。（不要被这样一个事实误导，即您可以通过指定tspan两个以上的元素来获得固定步长的输出——Matlab 的所有求解器都是可变步长求解器，并根据错误标准选择它们的真实步长。）

更好的选择是分段集成系统并将每次运行的结果输出附加在一起：

% t = 0 to t = 10, pass parameter a = 0 to add to ODEs
a = 0;
tspan = [0 10];
[T,Y] = ode45(@(t,y)myfun(t,y,a),tspan,y_0);

% t = 10 to t = 20, pass parameter a = 10 to add to ODEs
a = 10;
[t,y] = ode45(@(t,y)myfun(t,y,a),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

% t = 20 to t = 30, pass parameter a = 20 to add to ODEs
a = 20;
[t,y] = ode45(@(t,y)myfun(t,y,a),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

Matlab 编辑器可能会抱怨数组没有被预分配T和Y/或增长，但在这种情况下没关系，因为它们只增长了几次。或者，如果你想要固定的输出步长，你可以这样做：

dt = 0.01;
T = 0:dt:30;
Y = zeros(length(T),length(y_0));

% t = 0 to t = 10, pass parameter a = 0 to add to ODEs
a = 0;
[~,Y(1:10/dt+1,:)] = ode45(@(t,y)myfun(t,y,a),T(1:10/dt+1),y_0);

% t = 10 to t = 20, pass parameter a = 10 to add to ODEs
a = 10;
[~,Y(10/dt+1:20/dt+1,:)] = ode45(@(t,y)myfun(t,y,a),T(10/dt+1:20/dt+1),Y(10/dt+1,:));

% t = 20 to t = 30, pass parameter a = 20 to add to ODEs
a = 20;
[~,Y(20/dt+1:end,:)] = ode45(@(t,y)myfun(t,y,a),T(20/dt+1:end),Y(20/dt+1,:));

如果需要，可以轻松地将上述两个代码块转换为更紧凑的for循环。

在这两种情况下，您的 ODE 函数都会以这种方式myfun合并参数a：

function ydot = myfun(t,y,a)
    y(1) = ... % add a however you like
    ...
于 2013-06-28T17:48:08.823 回答
0

好的，就像 Simon McKenzie 说的那样，我们确实需要有关您的紧急问题的更多信息，但我想我可以提供帮助。根据您给我们的信息，我假设您有一个myfun传递给类似的函数ode45

y_0 = [0,0,0];
% here Tfinal is the time in seconds that you want to simulate to
% and specify the tspan so that you will get a solution at each whole
% number of seconds, I believe
[t,y] = ode45(@myfun,[0:0.1:Tfinal],y_0);

你在某个地方定义了你的函数，在这里我会调用它myfun

function dy = myfun(t,y)
  % here let's check to see if it's time to add your offsets
  % still, you may want to put a little fudge factor for the time, but
  % you'll have to experiment, I'll set it up though
  EPS = 1e-12;
  if( t < 10 + EPS || t > 10 - EPS )
    y(1) = y(1) + 10;
  .
  .
  .
  .
  % this only makes sense if you then use y(1) in the compuatation

否则，只需将偏移量添加到返回的解向量中，即

idx10 = find( t == 10 ); % should be there since it's in the tspan
y(idx10:end) = y(idx10:end) + 10; % I guess you add it from that point on?
于 2013-06-28T13:48:08.693 回答