How solve a system of ordinary differntial equation with time-dependent parameters

Question

How solve a system of ordinary differential equation ..an initial value problem ....with parameters dependent on time or independent variable? say the equation I have

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

where a(t) is a vector and b(t) =c*a(t); where the value of a and b are changing with time not in monotone way and each time step. I tried to solve this using this post....but when I applied the same principle ...I got the error message

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

Can someone please help me out?

Answer 1

Please read until the end to see whether the first part or second part of the answer is relevant to you:

Part 1: First create an .m file with a function that describe your calculation and functions that will give a and b. For example: create a file called fun_name.m that will contain the following code:

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

Then use a second file with the following code:

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

Of course for the fun_name.m case I wrote you need not use sub functions for a(t) and b(t), you can just use the explicit functional form in Dy if that is possible (like cos(t) etc).

Part 2: If a(t) , b(t) are just vectors of numbers you happen to have that cannot be expressed as a function of t (as in part 1), then you'll need to have also a time vector for which each of them happens, this can be of course the same time you'll use for the ODE, but it need not be, as long as an interpolation will work. I'll treat the general case, when they have different time spans or resolutions. Then you can do something of the following, create the fun_name.m file:

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)] ;

In order to use it, see the following script:

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