3

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?

dinkelk
  • 2,676
  • 3
  • 26
  • 46
prashanta_himalay
  • 149
  • 2
  • 2
  • 9
  • 2
    I assume you are using `ode45`? can you show the MATLAB function you are integrating, and how you are calling the integrator? – dinkelk Jan 15 '13 at 18:40

1 Answers1

4

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');
bla
  • 25,846
  • 10
  • 70
  • 101
  • So far so good and thanks a lot for your detailed answer.But another complication arises .... and That is actually the problem in my code.In my code,parameters will be come zero after certain point of time i.e say from 0=0 && t<=tau)a=interp1(taun,B1t,t) else a=0 I get absolutely wrong answer as the function is taking the t vector not individual t and corresponding a. – prashanta_himalay Jan 16 '13 at 13:03
  • that's a different problem, I would use the conditions in the code for the ODE itself, similar to what is discussed here: http://stackoverflow.com/questions/14173377/matlab-tricky-ode-system-with-boolean – bla Jan 16 '13 at 16:44
  • @bla: the answer you linked to is a very poor one. Introducing such conditions in the ODE function is likely to make the system stiff, and possibly result in completely incorrect results. At a minimum `ode15s` should be used, but even that can fail in certain cases. The proper way to handle such cases is via event detection. See my answers to [1](http://stackoverflow.com/q/17356037/2278029), [2](http://scicomp.stackexchange.com/q/16234/4474), and [3](http://stackoverflow.com/q/28432963/2278029). FYI, I got here from [this recent post](http://stackoverflow.com/q/32340707/2278029). – horchler Sep 01 '15 at 22:07
  • Well, it's not clear whether you mean part 1 or part 2 of the answer. Also, if you are right, please answer this question in the way you see fit. We can compare a simple case and see if my answer is satisfactory. I'll be happy to see and learn of other ways – bla Sep 02 '15 at 00:07