以您建议的方式(以及@macduff所说明的方式)在您的ODE中插入较大的不连续性可能会导致较低的精度和较长的计算时间(特别是使用ode45 - ode15s可能是更好的选择，或者至少确保您的绝对和相对公差是合适的)。您已经有效地生成了一个非常stiff system。如果您想要在特定时间开始的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
    ...

票数 2

Stack Overflow用户

发布于 2013-06-28 21:48:09

好的，就像西蒙McKenzie说的那样，我们真的需要更多关于你的紧急问题的信息，但我想我可以提供帮助。根据您提供的内容，我假设您有一个传递给ode45之类的函数的myfun

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?

票数 0

页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持

原文链接：

https://stackoverflow.com/questions/17356037

复制

相似问题

问Matlab ode解算器:更改状态和指定时间
EN

回答 2

Stack Overflow用户

Stack Overflow用户

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问Matlab ode解算器:更改状态和指定时间EN

回答 2

Stack Overflow用户

Stack Overflow用户

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问Matlab ode解算器:更改状态和指定时间
EN