matlab - 如何求解具有 transient 参数的常微分方程组

Question

如何求解参数依赖于时间或自变量的常微分方程组......初始值问题......？ 说出我的等式

 Dy(1)/dt=a(t)*y(1)+b(t)*y(2);
 Dy(2)/dt=-a(t)*y(3)+b(t)*y(1);
 Dy(3)/dt=a(t)*y(2);

其中 a(t) 是向量，b(t) =c*a(t)；其中 a 和 b 的值随时间变化，而不是单调的方式和每个时间步长。 我尝试使用 this post 来解决这个问题....但是当我应用相同的原理时...我收到了错误消息

"Error using griddedInterpolant The point coordinates are not sequenced in strict monotonic order."

有人可以帮帮我吗？

Answer 1

请阅读到最后，看看答案的第一部分或第二部分是否与您相关:

第 1 部分: 首先创建一个 .m 文件，其中包含描述您的计算的函数和将给出 a 和 b 的函数。例如:创建一个名为 fun_name.m 的文件，其中将包含以下代码:

function Dy = fun_name(t,y)

Dy=[ a(t)*y(1)+b(t)*y(2); ...
    -a(t)*y(3)+b(t)*y(1); ...
     a(t)*y(2)] ;
end

    function fa=a(t);
        fa=cos(t); % or place whatever you want to place for a(t)..
    end

    function fb=b(t);
        fb=sin(t);  % or place whatever you want to place for b(t)..
    end

然后使用带有以下代码的第二个文件:

t_values=linspace(0,10,101); % the time vector you want to use, or use tspan type vector, [0 10]
initial_cond=[1 ; 0 ; 0];  
[tv,Yv]=ode45('fun_name',t_values,initial_cond);
plot(tv,Yv(:,1),'+',tv,Yv(:,2),'x',tv,Yv(:,3),'o');
legend('y1','y2','y3');

当然，对于我写的 fun_name.m 案例，您不需要为 a(t) 和 b(t) 使用子函数，如果可能的话，您可以在 Dy 中使用显式函数形式(例如 cos(t) 等)。

第 2 部分: 如果 a(t) 和 b(t) 只是您碰巧拥有的无法表达的数字向量作为 t 的函数(如第 1 部分)，那么您还需要有一个时间向量，它们中的每一个都会发生，这当然可以与您将用于ODE，但不需要，只要插值有效即可。当他们有不同的时间跨度或解决方案时，我将处理一般情况。然后您可以执行以下操作，创建 fun_name.m 文件:

function Dy = fun_name(t, y, at, a, bt, b) 

a = interp1(at, a, t); % Interpolate the data set (at, a) at times t
b = interp1(at, b, t); % Interpolate the data set (bt, b) at times t

Dy=[ a*y(1)+b*y(2); ...
    -a*y(3)+b*y(1); ...
     a*y(2)] ;

为了使用它，请看下面的脚本:

%generate bogus `a` ad `b` function vectors with different time vectors `at` and `bt`

at= linspace(-1, 11, 74); % Generate t for a in a generic case where their time span and sampling can be different
bt= linspace(-3, 33, 122); % Generate t for b    
a=rand(numel(at,1));
b=rand(numel(bt,1));
% or use those you have, but you also need to pass their time info...

t_values=linspace(0,10,101); % the time vector you want to use 
initial_cond=[1 ; 0 ; 0];  
[tv,Yv]= ode45(@(t,y) fun_name(t, y, at, a, bt, b), t_values, initial_cond); % 
plot(tv,Yv(:,1),'+',tv,Yv(:,2),'x',tv,Yv(:,3),'o');
legend('y1','y2','y3');