Solving System of ODEs using for cycles to reproduce summing

Dear all,

Thank you first all for your attention reading this post.

I am a kind of new in Matlab so it might be that the answer to my question is really easy but I don't see the solution for the moment.

I want to solve a system of coupled differential equations that in mathematical form are presented below.

I am creating a function for this system that looks like

function dXdz = coupled_ODE(z,X,y,n,alpha,beta,J,Pf,uf)

m = length(X);                  % Defines the number of elements of the independent variable
dXdz = zeros(m,1);              % Initializes the column vector of the elements to the left of the ODE 

for i = 1:n
    dXdz(1) = dXdz(1) - alpha*J(i)*(Pf*X(i+1) - X(m)*y(i));
end

for i = 1:n
    dXdz(i+1) = (-1/X(1))*(X(i+1)*dXdz(1) + alpha*J(i)*(Pf*X(i+1) - X(m)*y(i)));
end

dXdz(m) = beta*(uf - X(1))/X(m);

where the parameters from y to uf are given as inputs, since after each iteration in the main code those are updated. Notice also that alpha and beta are constants.

This function is called in the main body as

[Z, Xout] = ode45(@coupled_ODE,zspan,Xin,[],yi,num_species,alpha,beta,J,Pf,uf);

I am getting problems in the results I obtain. They are not at all what I expect. The output of this ode45 is always the same at each iteration, indicating that no numerical integration of the ODEs is being performed at all.

I don't see where the problem can be. I might suspect that is coming from this for cycles I introduced to reproduce the summings, but I am not sure. Any idea?

Hopefully, you can see the error (and I would be not surprised if it is an amateur mistake).

Thank you all! Greetings

[rantunes]

I just want to add that the input vector X is a sex-elements column vector, with the elements

X = [u xi p]';

where xi is a vector of 4 elements. This explains why I used "X(i+1)", for instance.

Answer 1

[Torsten]

Change the call to ODE45 to

[Z,Xout]=ode45(@(z,X)coupled_ODE(z,X,yi,num_species,alpha,beta,J,Pf,uf),zspan,Xin);

Best wishes

Torsten.

[rantunes]

I tried out your suggestion but the output is exactly the same. I add here some more code related with the integration interval:

      Xin = [u xi p]';
      zspan = [incr incr+h];
     [Z,Xout]=ode45(@(z,X)coupled_ODE(z,X,yi,num_species,alpha,beta,J,Pf,uf),zspan,Xin);

      u = Xout(end,1);
      for i = 1:num_species
          xi(i) = Xout(end,i+1);
      end
      p = Xout(end,end);

zspan is an interval that is increasing after each while loop iteration (in the end of the cycle, incr = incr + h) until a certain length L.

The values obtained from the ode are stored as new u, xi and p, that are given as input for the next iteration, which influence the values of yi and so on.

I don't really understand where might be the flaw.

[Torsten]

How do you store the solution for the plot over [0:L] ?

Why do you call ode45 in a while loop with increment h ? Why don't you let ODE45 integrate over [0:L] in one call ?

Best wishes

Torsten.

[rantunes]

The point is that I can't integrate from 0 to L at once since at each increment h, the parameters that influence the system of ODE's also changes (namely yi, and beta). Am I clear? So I guess I have to solve the ODE in a cycle or not?

Best. Rodrigo

[Torsten]

You can't update the parameters in the "coupled_ODE" function, depending on the solution variables u,xi and p ?

Best wishes

Torsten.

[rantunes]

But this wouldn't bring me any advantage since I need the output from the ODE to re-calculate those parameters that are afterwards used to calculate yi, beta and so on... So I need a cycle that uses the values from Xout, I guess. So the system of ODEs has to be solved iteratively.

Best. Rodrigo

[Torsten]

ODE45 usually works with adaptive stepping in z-direction. If you force it to take a step of at most h, you will waste simulation time and loose precision of the solution.

This might explain why you get constant output, namely if your "h" is chosen too small. Another reason might be that you don't store the solution correctly for the plot over [0:L].

Best wishes

Torsten.

[rantunes]

Yes of course. But even when I increase h, the same output is given. This is how the while cycle looks like:

incr = 0;
iter = 1;

while incr < b
    if incr == a                                                                 % If incr == a, the variables assume the value of the initial parameters                
        u = uf;
        yi = yini;
        xi = xfi;
    end

      yi = local_permeate_concentration(yi,J,xi,p,Pf);
      mu_mix = viscosity(yi,sigma,omega,M,T);                                         % Computation of viscosity of the mixture     
      beta = 128*R*T*mu_mix/(pi*((Di*1e-6)^4)*Nfiber);                                % Determination of the parameters needed to calculate the pressure profile (input coupled_ODE)
      summ = summing(J,xi,yi,Pf,p);                                                   % Calculation of the sum J(i)*Pf(x(i) - p*x(i))
      y_bulk = bulk_permeate_concentration(J,Pf,p,xi,yi,summ);                        % The bulk concentration with is computed

      Xin = [u xi p]';
  %     disp(num2str(Xin));
      zspan = [incr incr+h];
      [Z, Xout] = ode45(@coupled_ODE,zspan,Xin,[],yi,num_species,alpha,beta,J,Pf,uf);      % u, p and xi are computed using the initial values, with yi as input 
     % [Z,Xout]=ode45(@(z,X)coupled_ODE(z,X,yi,num_species,alpha,beta,J,Pf,uf),zspan,Xin);
      disp(num2str(Xout));

      u = Xout(end,1);
      for i = 1:num_species
          xi(i) = Xout(end,i+1);
      end
      p = Xout(end,end);
      %break;

      %y_bulk = bulk_permeate_concentration(J,Pf,p,xi,yi,summ);                        

      yi_save(iter,:) = yi;
      y_bulk_save(iter,:) = y_bulk;
      xi_save(iter,:) = xi;
      p_save(iter,1) = p;
      u_save(iter,1) = u;

      iter = iter + 1;
      incr = incr + h;
  end

I dont see where I am storing the solution in a not proper way... Until the ode, the functions are properly working. It is really how the ode45 is implemented, somehow.

Greets. Rodrigo

[Torsten]

Are you sure h is not equal to 0 ? I don't see where it is set.

If h were 0, u, yi and xi were always set to their initial values in the while-loop.

Please check Xout(end,:) of iteration j and the initial values Xin with which you enter ODE45 in iteration j+1. Are they identical ?

Best wishes

Torsten.

[rantunes]

h is equal to 0.01 (is defined before those instructions on the top of the last code). iter is increasing in steps of h as supposed. I can check it.

The Xout(end,:) after (iteration j) and before (iteration j+1) ode45 is also the same... so the values are being well manipulated, I would say.

Best. Rodrigo

[Torsten]

Then I don't see anything obviously wrong in your code. I'd check the dXdz values in coupled_ODE ; they should always be approximately zero during the integration if the solution does not change with z.

Changing the initial stepsize to a value at least smaller than h and strengthing the tolerances might also be an option.

I strongly recommend putting all the updates of the parameters done in the while-loop in coupled_ODE instead and integrate from z=0 to z=L in one step.

Best wishes

Torsten.

[rantunes]

I tried to do what you suggested about the stepsize and it really didn't work. It is true that the dXdz values are small, but even multiplying them by very big numbers the things don't really change. I guess it has to do how I implement the system, the system ODE. I have the feeling that using this for cycle to write those equations is not the optimal one...

The algorithm is in such a way that the parameters like yi at each step are used to integrate the ODE. So it is not really possible to integrate over 0:L outside the loop.

I really dont know what to do...

Best. Rodrigo

[Torsten]

Well, I don't see why the subseqent calls can not be done within coupled_ODE

yi = local_permeate_concentration(yi,J,xi,p,Pf);
        mu_mix = viscosity(yi,sigma,omega,M,T);                                         % Computation of viscosity of the mixture     
        beta = 128*R*T*mu_mix/(pi*((Di*1e-6)^4)*Nfiber);                                % Determination of the parameters needed to calculate the pressure profile (input coupled_ODE)
        summ = summing(J,xi,yi,Pf,p);                                                   % Calculation of the sum J(i)*Pf(x(i) - p*x(i))
        y_bulk = bulk_permeate_concentration(J,Pf,p,xi,yi,summ);

which would result in an integration over [0:L].

But if you really have to integrate in constant steps of deltaz=0.01, there is no advantage using ODE45 instead of including the integration directly in the while-loop.

Use an explicit Euler-method with stepsize deltaz=0.01 for u, xi and p and see what happens.

Best wishes

Torsten.

[rantunes]

I could put those calls inside the couple_ODE, but the point is that at iteration j we obtain u, xi and p, that are used in iteration j+1 to obtain yi,etc., and so on. So the things have to run sequentially and iteratively. I can't just separate.

I will try out the explicit Euler-method and see what I can get.

Best. Rodrigo

[Torsten]

Don't think in iterations, think in states.

If the parameters yi,.... can be calculated from the states of u, xi and p (which are given variables in coupled_ODE), you can calculate dXdz.

Best wishes

Torsten.

[rantunes]

Finally I figured out what was the problem. It was a value that I was considering wrong in terms of units that were providing very low dXdz and thus no "visible" integration.

The implementation above was always correct.

Thanks. Rodrigo