3

如何解决常微分方程系统..初始值问题....参数取决于时间或自变量？说出我的方程式

 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 的值随时间变化，而不是以单调方式和每个时间步长变化。我试图用这篇文章解决这个问题......但是当我应用相同的原理时......我收到了错误消息

“使用 griddedInterpolant 出错点坐标未按严格单调顺序排序。”

有人可以帮我吗？

4

1 回答 1

4

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

第 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');

当然，对于我写的情况，你不需要为andfun_name.m使用子函数，如果可能的话，你可以使用显式函数形式（比如等）。a(t)b(t)Dycos(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');
于 2013-01-15T22:58:27.843 回答