How to code an ODE System with paramters depending on variables

Hey guys!

I'm trying to simulate a system of differential equations. For that I made the function:

function dYdt = System(t,Y)

    i_L = Y(1);
    u_cp = Y(2);
    u_cs = Y(3);
    u_A = Y(4);
        
    C_A = 28e-9;
    C_p = 251e-12;
    C_s = 390e-12;
    R = 400;
    L = 370e-6;
    I_Load = 4;
    
    di_Ldt = (u_e-u_cs-u_cp)/L;
    du_cpdt = (1/C_p)*(i_L-a*((i_L*C_A+C_p*(I_Load+(u_A/R))*sign(u_cp))/(C_p+C_A)));
    du_csdt = i_L/C_s;
    du_Adt = abs(i_L)*a*(C_A/(C_A+C_p))-(I_Load-(u_A/R))*(1-a*(C_A/(C_A+C_p)));
    
    dYdt = [di_Ldt; du_cpdt;  du_csdt; du_Adt];
    
end

So far so good, but as you might have noticed the variables a and u_e are not declared anywhere.

The problem is, that 'a' depends on 'i_L', 'u_cp' and 'u_cs' and 'u_e' depends on 'i_L'. For that i made functions for each variable.

function aa = a_Function(i_L, u_cp, u_cs)

        if i_L == 0
            aa = 0;
        end

        if abs(u_cp) == u_cs
            aa =1;
        end    
    end

and

    function u_e = u_e_Function(i_L)

        if i_L >= 0
            u_e = 400;
        else
            u_e = -400;
        end
    end

I am not sure where to put those functions. Do I put those inside the System- function or I declare them in the main file, where I execute the Simulation. Because to get the variables a and u_e I need first to put the ODE- System through a Solver (for example ode45), but for that die ODE's need a and u_e.

Does anyone have any suggestions?

Answer 1

[Alan Stevens]

You can include the body of each "extra" function directly in the main function:

function dYdt = System(t,Y)

    i_L = Y(1);
    u_cp = Y(2);
    u_cs = Y(3);
    u_A = Y(4);
        
    C_A = 28e-9;
    C_p = 251e-12;
    C_s = 390e-12;
    R = 400;
    L = 370e-6;
    I_Load = 4;
    
    if i_L == 0
        a = 0;
    end
    if abs(u_cp) == u_cs
        a =1;
    end  
    if i_L >= 0
         u_e = 400;
    else
         u_e = -400;
    end
    
    di_Ldt = (u_e-u_cs-u_cp)/L;
    du_cpdt = (1/C_p)*(i_L-a*((i_L*C_A+C_p*(I_Load+(u_A/R))*sign(u_cp))/(C_p+C_A)));
    du_csdt = i_L/C_s;
    du_Adt = abs(i_L)*a*(C_A/(C_A+C_p))-(I_Load-(u_A/R))*(1-a*(C_A/(C_A+C_p)));
    
    dYdt = [di_Ldt; du_cpdt;  du_csdt; du_Adt];
    
end

[Aleks Boehler]

Thank you for your answer. I did that earlier and the code works but the result is wrong and I think the problem ist that when ODE45 solves the equation it sets the values of a and u_e ones but doesnt update it each time the other variables change their values.

[Alan Stevens]

As long as the values in Y update each time, then a and u_e will update also.

To check, display the values of Y, a and u_e within the function each time it's called (i.e. remove the semicolons from i_L etc). Do they change?

[Aleks Boehler]

Ok then, that wasn't clear to me. Thank you very much.

Then the problem must be something outside of the coding.

I would like to check if the values change but now for some reason I get the message 'Unrecognized function or variable 'a'.' Even though the statement inside the if command is clearly there and when I mark 'a' in the equation with the courser the other 'a's in the if- command show up blue too. So I clearly didnt spell this one letter wrong :D

[Aleks Boehler]

Plus, everytime I run the System function I get the message that I have not enough function inputs. But when the solver uses the function it works... well it did up to the point where it suddenly decided not to recognice 'a' anymore.

[Alan Stevens]

Difficult for me to help further unless you upload all the code that calls the System function.

[Aleks Boehler]

tRange = [0, 20e-6];
Y0 = [4;0;0;0];

[tSol, YSol] = ode45(@System, tRange, Y0);

I_L = YSol(:,1);
U_cp = YSol(:,2);
U_cs = YSol(:,3);
U_A = YSol(:,4);



figure(1)
plot(tSol,I_L,'ko')
hold on 
plot(tSol,U_cp)
plot (tSol, U_cs)
hold off

figure(2)
plot(tSol, U_A)

and here comes the function ones again:

function dYdt = System(t,Y)

    i_L = Y(1);
    u_cp = Y(2);
    u_cs = Y(3);
    u_A = Y(4);
   
    
    if i_L == 0
         a = 0;
    end
       
    if abs(u_cp) == u_cs
        a = 1;
    end
    
    if i_L >= 0
        u_e = 400;
    else
        u_e = -400;
    end 
    
    
    a
    u_e
    C_A = 28e-9;
    C_p = 251e-12;
    C_s = 390e-12;
    R = 400;
    L = 370e-6;
    I_Load = 4;
    
    di_Ldt = (u_e-u_cs-u_cp)/L;
    
    du_cpdt = (1/C_p)*(i_L-a*((i_L*C_A+C_p*(I_Load+(u_A/R))*sign(u_cp))/(C_p+C_A)));
    du_csdt = i_L/C_s;
    du_Adt = abs(i_L)*a*(C_A/(C_A+C_p))-(I_Load-(u_A /R))*(1-a*(C_A/(C_A+C_p)));
    
    dYdt = [di_Ldt; du_cpdt;  du_csdt; du_Adt];
    
end

[Alan Stevens]

Ok. You need to specify default values for a and u_e. Your if..else..end options didn't necessarily cover all possibilities! Like so:

function dYdt = System(~,Y)

    i_L = Y(1);
    u_cp = Y(2);
    u_cs = Y(3);
    u_A = Y(4);
   
    a = 0;  % Specify a default for a 
%     if i_L == 0
%          a = 0;
%     end
    if abs(u_cp) == u_cs
        a = 1;
    end
    
    u_e = 400;  % Specify a default for u_e
    if i_L < 0
       u_e = -400; 
    end
%     if i_L >= 0
%         u_e = 400;
%     else
%         u_e = -400;
%     end 
    
    
    
    C_A = 28e-9;
    C_p = 251e-12;
    C_s = 390e-12;
    R = 400;
    L = 370e-6;
    I_Load = 4;
    
    di_Ldt = (u_e-u_cs-u_cp)/L;
    
    du_cpdt = (1/C_p)*(i_L-a*((i_L*C_A+C_p*(I_Load+(u_A/R))*sign(u_cp))/(C_p+C_A)));
    du_csdt = i_L/C_s;
    du_Adt = abs(i_L)*a*(C_A/(C_A+C_p))-(I_Load-(u_A /R))*(1-a*(C_A/(C_A+C_p)));
    
    dYdt = [di_Ldt; du_cpdt;  du_csdt; du_Adt];
    
end

It now works ok. Can't say if the output is sensible - that's for you to decide!

[Aleks Boehler]

With the default variable it works for me too. The rest I have to figure out by my self. Thank you very much once again. Did you also got the message 'Not enough Input paramaters' when you try to run the function?

[Alan Stevens]

One other thing I've noticed:

Your u_cp and u_cs are both floating point values. It's best not to do equality tests for these because of small rounding (or other) errors. Best to replace the absolute value of the difference with a small tolerance.

i.e. replace

if abs(u_cp) == u_cs
        a = 1;
    end

with

if abs(u_cp-u_cs) < tol
        a = 1;
    end

where you replace the word 'tol' with a small number (e.g. 10^-3 or 10^-6 or ... whatever is appropriate for your system).

[Alan Stevens]

I didn't get the 'not enough input parameters' message.