Pass out a dependent variable from ode15s without taking the derivative of it

I have search for hours but i can't seem to figure out how it is possible to exstract a variable depending on x (x is derrived using the ode15s)

My solution to a system is (simplified)

function [y,t,rr] = simulation(t,initial)
    options = odeset('RelTol', 1e-5, 'AbsTol', 1e-7);   
    [y,t,rr] = ode15s(@kinetics, t, initial, options); 
end

The kinetics function is as follows (simplified)

function [dxdt,x,rr] = kinetics(t,x,init)
%code desribing i different phenomenas all depending on "x" e.g:
Ph1 = numaxG * x(1,1) * x(8,1) / (KSPG + x(1,1) + (x(1,1)^2)/KiPG);
%constants like numaxG has been defined..

rr(1,1)  = Ph1 * Ph5 * Ph7 * Ph10 * Ph16 * Ph19a;                      
rr(2,1)  = Ph2 * Ph6 * Ph8 * Ph11 * Ph17 * Ph19b * Ph21;          
rr(3,1)  = Ph3;                                                       
rr(4,1)  = Ph12 *  Ph9;                                                
rr(5,1)  = Ph14;                                              
rr(6,1)  = 0;                                                          
rr(7,1)  = 0;

% some other code describing dxdt
end

right now rr is also integrated using the solver, but i would like it to get as an ouput without it being solved. For me this is difficult as it is dependent on x which is derived.

[darova]

I think you are using ode15s in a wrong way. Description from help page

[Andreas Skovhøj]

im not sure what you mean.. i just use it this way. What should "te" be useful for?

[Torsten]

You used three output parameters in your code, not two.

[Andreas Skovhøj]

no but that is my problem.. i know that i should only use two but i added "rr" because i need it as a output from the function.

Answer 1

[darova]

You should calculate rr variable after solving ode45

An example

function main
[t,u] = ode45(@f,[0.1 3],[1 1]);
u1 = u(:,1);        % extract u1
u2 = u(:,2);        % extract u2
rr1 = u1 - u2 ...
end

function du = f(t,u)
du(1,1) = u(2);
du(2,1) = sin(u(1))/t;
end

[Andreas Skovhøj]

Dear darova,

Thank you so much for your help.

As you also write in your example: "u" is the result from solving the function "f" using the ode45 but that is not the value I want to exstract. The value i want to exstract is (in my case) the x vector which is changing in every step using the ODE45. So how do exstract x for every step and use this value to calculate my rr vector.

As an example: Right now if i type "rr" in my "kinetics" function which i solve using ode45, i can see the rr vector for every step in the ode45 when i run the code. And THAT is the rr vector which i would like to store so i can plot it later.

Hope it makes some kind of sense

Regards

Andreas

[darova]

• So how do exstract x for every step and use this value to calculate my rr vector.

[t,u] = ode45(@f,[0.1 3],[1 1]);
u(:,1)      % x1
u(:,2)      % x2

[Andreas Skovhøj]

hmm i dont know if we are on the same page (im not an experienced MATLAB user). can you maybe give another example where a vector (depending on x) inside your function @f is ignored by the solver, but just calculated using the given x value of each step)

[darova]

Can you attach original equations?

[darova]

Assume you have equations (t - time, x=x(t) - unknown function)

solution for this problem

x0 = 1;                 % initial condition
dx = @(t,x) x^2-t;
[t,x] = ode45(dx,[0 2],x0);
r = sin(x)

[Andreas Skovhøj]

Here are my two files i think you will need to understand by problem

[darova]

Here is an idea

Your current kinetics.m

function [dxdt,rr,x] = kinetics(t,x,p,init,t_load,a,mode)
    % long code
    % computations with 'x' and other parameters
    rr(1:7,1) = ... % calculations
    dxdt = ...      % derivatives
end

Write an additional function

function par = mycalc(p,x)
    % long code
    % computations with 'x' and other parameters
    % ... return all needed results for further calculation as struct
    par.Ph1 = Ph1;
    par.Ph2 = Ph2;
end
function dxdt = kinetics(t,x,p,init,t_load,a,mode)
    par = mycalc(p,x)
    rr(1:7,1) = ... % calculations
    dxdt = ...      % derivatives
end

And after solving (having x and all needed parameters) you can just

par = mycalc(p,x);
rr = par.rr;

I don't like to use structures since each parametr should be written (you have a lot of them). But it's the only idea i have for your case

I saw this line at the end

dxdt = rr'*st + vol';

Once you solved your equations you have x. You can calculate dxdt (using diff possibly)

So it possible to calculate rr

rr = inv(st)*(dxdt-vol)