Solving a set of differential equations using ode45 with a time dependent parameter

Hello,

I'm trying to solve a set of differential equations using ode45 (code below). The parameter 'Flux' is time dependent. At this moment the ModelFotoconversie_ODE function creates an array with the values of 'Flux' at each time. How can I adapt my code in such way that the ModelFotoconversie_ODE function only calculates the parameter 'Flux' at the specific time t ? _________________________________________________________________________________________________________ function [t,C] = ODE(t, C_SP0,C_MC0,pars,Vtot,L,eps,Flux,Nu,RatioDark)

tspan = [0 t(end)];

y0 = [C_SP0 C_MC0];

[t, C] = ode45(@(t,C) ModelFotoconversie_ODE(t,C,pars,Vtot,L,eps,Flux,Nu,RatioDark), tspan, y0);

end

_________________________________________________________________________________________________________

function [ DC ] = ModelFotoconversie_ODE(t,C,pars, Vtot,L,eps,Flux,Nu,RatioDark)

k(1,1)=pars(1);

k(1,2)=k(1,1)*RatioDark;

QE(1,1)=pars(3);

QE(1,2)=pars(4);

N_Av= 6.022140857e23;

h = 6.626*10^(-34); %J.s

Abs_tot = [eps(1,1)*C(1)+eps(1,2)*C(2)' , eps(2,1)*C(1)+eps(2,2)*C(2)'];

Flux_Transm = Flux.*10.^(-Abs_tot.*L);

N_abs_tot = ((Flux- Flux_Transm)./(h*Nu))/1000;

N_Abs_UV = [eps(1,1)*C(1) eps(1,2)*C(2)].*(N_abs_tot/Abs_tot(1,1));

N_Abs_Geel = [eps(2,1)*C(1) eps(2,2)*C(2)].*(N_abs_tot/Abs_tot(1,2));

N_Abs= [N_Abs_UV; N_Abs_Geel] ;

N_Abs_SP= sum(N_Abs(:,1));

N_Abs_MC= sum(N_Abs(:,2));

dC1= k(1,2)*C(2) - k(1,1)*C(1)+(QE(1,2)*N_Abs_MC - QE(1,1)*N_Abs_SP)/(Vtot*N_Av);

dC2= k(1,1)*C(1) - k(1,2)*C(2)+(QE(1,1)*N_Abs_SP - QE(1,2)*N_Abs_MC)/(Vtot*N_Av);

DC = [dC1; dC2]; end

Answer 1

[Star Strider]

In the ode45 (link) documentation, see the section on ODE with Time-Dependent Terms (link).

[Bram Callewaert]

Hello,

First of all thanks for the quick answer. I already read the documentation before i asked the question, but i don't understand how i can use this in my case.

I made a linspace in the ODE function:

ft = linspace(0,340647,340648);

But i dont see how i have to make f

[Star Strider]

I don’t know what ‘Flux’ is, although I assume it’s a vector. The documentation uses the interp1 function to interpolate a value for the time-dependent variable at a particular time. The time is an argument to the ODE function, so use that as the query value in the interp1 call to return the value of ‘Flux’ at ‘t’.

Note that for this, ‘Flux’ is defined inside the ODE function, not passed as a parameter, and the interp1 call is simply used to provide an appropriate value for it.

[Bram Callewaert]

'Flux' is a 19x2 vector.

if i adapt the code like this:

function [t,C] = ODE(t, C_SP0,C_MC0,pars,Vtot,L,eps,Flux,Nu,RatioDark)

tspan = [0 t(end)];

y0 = [C_SP0 C_MC0];

ft = linspace(0,340647,340648);

[tijd, C] = ode45(@(t,C) ModelFotoconversie_ODE(t,C,pars,Vtot,L,eps,Flux, ft,Nu,RatioDark), tspan, y0);

end

______________________________________________

function [ DC ] = ModelFotoconversie_ODE(t,C,pars, Vtot,L,eps,Flux,ft,Nu,RatioDark)

% k=[k_SP_to_MC k_MC_to_SP] % QE=[QE_SP_to_MC QE_MC_to_SP] % eps=[eps_SP_UV eps_MC_UV;eps_SP_VIS eps_MC_VIS] % % C=[C_SP C_MC]

k(1,1)=pars(1); k(1,2)=k(1,1)*RatioDark; QE(1,1)=pars(3); QE(1,2)=pars(4);

N_Av= 6.022140857e23; h = 6.626*10^(-34); %J.s

f = interp1(ft,Flux,tijd);

Abs_tot = [eps(1,1)*C(1)+eps(1,2)*C(2)' , eps(2,1)*C(1)+eps(2,2)*C(2)'];

Flux_Transm = f.*10.^(-Abs_tot.*L);

N_abs_tot = ((f - Flux_Transm)./(h*Nu))/1000;

N_Abs_UV = [eps(1,1)*C(1) eps(1,2)*C(2)].*(N_abs_tot/Abs_tot(1,1));

N_Abs_Geel = [eps(2,1)*C(1) eps(2,2)*C(2)].*(N_abs_tot/Abs_tot(1,2));

N_Abs= [N_Abs_UV; N_Abs_Geel];

N_Abs_SP= sum(N_Abs(:,1)); N_Abs_MC= sum(N_Abs(:,2));

dC1= k(1,2)*C(2) - k(1,1)*C(1)+(QE(1,2)*N_Abs_MC - QE(1,1)*N_Abs_SP)/(Vtot*N_Av);

dC2= k(1,1)*C(1) - k(1,2)*C(2)+(QE(1,1)*N_Abs_SP - QE(1,2)*N_Abs_MC)/(Vtot*N_Av);

DC = [dC1; dC2];

end

_______________________________________________

I get the following error message:

Error using griddedInterpolant The grid vectors do not define a grid of points that match the given values.

what am i doing wrong ?

[Star Strider]

I can’t understand your code, and without formatting (use the [{}Code] button), it’s difficult to read.

What information do you want to get out of ‘Flux’ at time ‘t’? I don’t see an assignment for ‘tijd’.

That seems to be the issue, so let’s concentrate on that just now.

[Bram Callewaert]

Im sorry i replaced all the 'tijd' in the previous example but forgot that one. This is what is should look like: ( I would like to have the Value of Flux at time 'tijd')

The vectors 'tijd' and Flux have the same length. Flux contains the flux value at each time 'tijd'.

function [tijd,C] = ODE(tijd, C_SP0,C_MC0,pars,Vtot,L,eps,Flux_metingen,Nu,RatioDark)

    tspan = [0 tijd(end)];
    y0 = [C_SP0 C_MC0];
    ft = linspace(0,340647,340648);

    function [ DC ] = ModelFotoconversie_ODE(tijd,C,pars, Vtot,L,eps,Flux_metingen,ft,Nu,RatioDark)

end

___________________________________________________________

function [ DC ] = ModelFotoconversie_ODE(tijd,C,pars, Vtot,L,eps,Flux_metingen,ft,Nu,RatioDark

    k(1,1)=pars(1);
    k(1,2)=k(1,1)*RatioDark;
    QE(1,1)=pars(3);
    QE(1,2)=pars(4);
    N_Av= 6.022140857e23;
    h = 6.626*10^(-34); %J.s

    f = interp1(ft,Flux_metingen,tijd);

    Abs_tot = [eps(1,1)*C(1)+eps(1,2)*C(2)' , eps(2,1)*C(1)+eps(2,2)*C(2)']; 

    Flux_Transm = f.*10.^(-Abs_tot.*L);
    N_abs_tot = ((f - Flux_Transm)./(h*Nu))/1000;
    N_Abs_UV = [eps(1,1)*C(1)  eps(1,2)*C(2)].*(N_abs_tot/Abs_tot(1,1)); 
    N_Abs_Geel = [eps(2,1)*C(1)  eps(2,2)*C(2)].*(N_abs_tot/Abs_tot(1,2)); 

    N_Abs= [N_Abs_UV; N_Abs_Geel] ;
    N_Abs_SP= sum(N_Abs(:,1));
    N_Abs_MC= sum(N_Abs(:,2));

    dC1= k(1,2)*C(2) - k(1,1)*C(1)+(QE(1,2)*N_Abs_MC - QE(1,1)*N_Abs_SP)/(Vtot*N_Av);
    dC2= k(1,1)*C(1) - k(1,2)*C(2)+(QE(1,1)*N_Abs_SP - QE(1,2)*N_Abs_MC)/(Vtot*N_Av);

    DC = [dC1; dC2];
  end

[Star Strider]

I’m more lost than I was before.

What is ‘Flux’ and what is the time vector that goes with it?

I was thinking of something like this:

Flux = [sin(2*pi*[0:18]/9); cos(2*pi*[0:18]/9)]';
tFlux = [0:18]';

t = 10.2;

Flux_at_t = interp1(tFlux, Flux, t)

Flux_at_t =
      0.71119      0.64757

[Bram Callewaert]

This is what 'Flux' looks like: Flux_metingen =

         0         0
         0         0
         0         0
         0   83.0000
         0         0
         0   83.0000
    4.7000         0
    4.7000         0
    4.7000         0
    4.7000         0
   60.0000         0
   60.0000         0
         0         0
         0    8.0000
         0    8.0000
         0   83.0000
         0   83.0000
         0         0
         0         0

The time vector tijd looks like :

1.0e+05 *

         0
    0.0093
    0.0167
    0.0218
    0.0335
    0.0399
    0.0424
    0.0451
    0.0477
    0.0508
    0.0556
    0.0606
    0.0648
    0.0677
    0.0708
    0.0736
    0.0789
    0.9461
    3.4065

So at time tijd = 1.0e+05 *0.0218 we have Flux [0 83], at time 1.0e+05 *0.0335 we have Flux [0 0]

this time vector tijd is also the vector wich is used to solve the differential equations.

What i am trying is to let the function ModelFotoconversie_ODE calculate the value of Flux at the same time from 'tijd' that it uses for that step of the differential equation.

I guess (correct me if i am wrong) i need something like:

tFlux = tijd; (in ODE)

and

f = interp1(tFlux,Flux_metingen,tijd); (in ModelFOtoconversie_ODE)

[Star Strider]

I would put ‘Flux’ and ‘tjid’ inside your ‘ModelFOtoconversie_ODE’ function, then do:

f = interp1(tjid, Flux_metingen, t, 'linear', 'extrap');

I added ‘'linear','extrap'’ to keep the interpolation from crashing if it inadvertently exceeds the minimum or maximum values of ‘tjid’.

[Bram Callewaert]

Now:

[tijd, C] = ode45(@(tijd,C) ModelFotoconversie_ODE(tijd,C,pars,Vtot,L,eps,Flux_metingen, tFlux,Nu,RatioDark), tspan, y0);

returns mostly NaN's

C =

   1.0e+04 *

    0.0000    0.0000
       NaN       NaN
       NaN       NaN
       NaN       NaN
    4.0154   -4.0154
       NaN       NaN
       NaN       NaN
       NaN       NaN
       NaN       NaN
       NaN       NaN
       NaN       NaN
       NaN       NaN
       ...       ...

[Star Strider]

The NaN values are the result of calculations producing 0/0, Inf/Inf or other values. One way to deal with that is to replace the 0 in ‘Flux’ with eps.

I cannot follow your code, but I notice that you have a variable named eps. Please rename your variable! Until you do, my fix for your NaN problem will not work.

[Bram Callewaert]

I renamed eps and chanced the zero's in Flux. But i am stil getting NaN.

I also tried to ad a + eps in N_Abs_UV and N_Abs_Geel to avoid a 0/0 in the division:

   N_Abs_UV = [epsilon(1,1)*C(1)  epsilon(1,2)*C(2)].*(N_abs_tot(1)/Abs_tot(1)+eps);

But it still isnt working

[Star Strider]

I can’t follow your code. It will be worthwhile for you to use the debug functions in the Editor to step through your code to see where the problems arise. That’s the only option I can think of at this point.