How can I stop an ODE overwriting variable values?

As part of my ODE I am calculating values that change with iteration, however they are not part of the ODE outputs. Could someones provide me with a way to stop these values being overwritten so I can create a graph of these values?

ODE

%^H+ concentration quadratic

a = 1;

b = C_Na_plus + K_D - C_C_neg - C_D_neg - C_Cl_neg;

c = C_Na_plus*K_D - K_w - K_D*C_A_Total- K_D*C_C_neg - K_D*C_D_neg - K_D*C_Cl_neg;

d = -K_w*K_D;

r = roots([a b c d]);

r = (r(isreal(r)&real(r)>=0));

C_H_plus = r

C_A_neg = ((K_D*C_A_Total)/(C_H_plus+K_D));

C_A = ((C_A_neg * C_H_plus)/K_D);

C_OH_neg =( K_w/C_H_plus)

K_D_Check = (C_A_neg*C_H_plus)/C_A;

C_A_Total = C_A + C_A_neg;

Hatta = sqrt(k1*C_A_neg*Dif)/k_L

% Define k2

k2 = 4.46*10^-6 + 6.018*C_OH_neg;

% Rate Equations & Mass transfer

r_NaC = k1 * C_A_neg * C_B;

r_NaD = k2 * C_B;

r_A2 = ((k3 * K_D * C_H_plus)/(K_D + C_H_plus)^2)*C_A^2;

C_BL = (1/H)*(1/(1+((k1*C_A_neg*El)/k_La)));

JB = k1*C_BL*C_A_neg*El;

% PH

pH = -log10(C_H_plus)

% ODES

dFdV(1) = -r_NaC-(2*r_A2);

dFdV(2) = -r_NaC -r_NaD + JB;

dFdV(3) = r_NaC;

dFdV(4) = r_NaD;

dFdV(5) = 0;

dFdV(6) = 0;

dFdV(7) = r_NaC +r_NaD;

dFdV(8) = r_NaC +r_NaD;

dFdV(9) = r_A2;

dFdV(10)= 0;

The aim is basically to also plot pH and Hatta with respect to V.

[Geoff Hayes]

JF - please format the above code that it is kept distinct from the rest of your question. Highlight the code portion and press the {} Code button.

As Hatta and pH are just initialized once, then it is unclear how it can be overwritten. Please include all of your code that may be relevant to the question. Is the above part of a body of a function that you call with different inputs or..? Where is C_Na_plus and the rest of the local variables defined/initialized? If you feel that there is too much code to paste above, then just attach the m file or files to your question using the paperclip button.

Answer 1

[Matt Tearle]

If I understand correctly, your ODE function takes in t and y (and possibly some other parameters), then calculates Hatta and pH from those, and then the derivatives. You would like the values of Hatta and pH that are calculated at each time step, right?

There's no immediate way to do that through ode45. The ODE function is actually called a bunch of times at a bunch of intermediate locations that aren't necessarily part of the output from ode45 (that's how RK methods work).

So basically you need to calculate Hatta and pH from the t and y that was returned by ode45. Some redundant effort, yes, but that's the way it is.

Probably the easiest way to do it is to make a function that returns Hatta and pH from t and y and call that in your ODE function. Then, once ode45 gives you back your solution t and y you can just call that function directly to get the corresponding values of Hatta and pH.

The tricky thing is that ode45 will be calling your ODE function with a scalar t and a column vector y. But the solution you get back (that you'll then pass to your intermediate function) is in the form of a column vector t and a matrix y. That can cause problems when you reference elements of y (such as y(3) vs y(:,3)).

Here's a simple example:

function dy = odefunc(t,y)
intval = y(2) - t; % I want this value at the end
dy = [y(1) + intval;y(1).*y(2)];

Script to run it:

[t,y] = ode45(@odefunc,[0 1],[0 1]);
plot(t,y(:,1),t,y(:,2))

Now, to get intval as a function of t, change the ODE function:

function dy = odefunc(t,y)
intval = odeintval(t,y);
dy = [y(1) + intval;y(1).*y(2)];

function intval = odeintval(t,y)
intval = y(2,:) - t;  % Note the use of matrix indexing (works for column vector or matrix)

And run the script:

[t,y] = ode45(@odefunc,[0 1],[0 1]);
plot(t,y(:,1),t,y(:,2))

intval = odeintval(t',y');  % Note the transposes (so second element of y is y(2,:) instead of y(:,2))
figure
plot(t,intval)