Using for loop to calculate ode15s multiple time, store end value of each calculation

I want to simulate the magma flow in the dike, and I want to calculate the ode function multiple time with changing parameters (changing magma discharge rate)

And get z and y(1), y(2), y(3), y(4), y(5) at end (event happened), by each Q value.

I tried using for loop to change "Q" each time, but it does not work, maybe matlab try to calculate Q as a whole array, not chaning parameter each time.

Could you give any suggestions, or advice?

sol = cell(101,1); %assigning end value
zspan = [-3200 1000]; %distance start from -3200
y0 = [79000000; 0.000001; 2.85; 79100000; 0.05; 0]; % initial value 
% constants
g = 9.8; a = 0.5; b = 200; Nb = 10^12; sigma = 0.05;
T = 1353; IFF = 10^(-13.5); kv = 4*pi/3; phim = 0.6; phimb = 0.8; 
rhol = 2500; Ro = 0.000001; B = 8.31;  Rw = 461;
So = (3/(4*pi))^(1/3)*(1/Nb)^(1/3); Mw = 0.018;

%% ode45 function
Opt1    = odeset('Events', @myEvent1);
Opt2    = odeset('Events', @myEvent3);

Opt = odeset('Events', @myEvent2);

Q = linspace(500,510,101); %Changing parameter Q for each loop
for k = 1:101
    sol(k) = ode15s(@(z,y) Tuto31(z,y,Q(k)),zspan,y0,Opt); %Assigning end value of fundction 
end
y = y(y==real(y));

function [value, isterminal, direction] = myEvent2(z, y); %Event function to end
phim = 0.45;
phieq = 0.0184*y(3).^3-0.043*y(3).^2-0.2238*y(3)+0.6691;
Nb = 10^12; Ro = 0.000001;
So = (3/(4*pi))^(1/3)*(1/Nb)^(1/3);
S = ((So.^3)-(Ro.^3)+(y(2).^3))^(1/3);
phib = (y(2)^3)/(S.^3);
if y(5)<0.4409
visxb = (1-(y(5)/phim)).^(-5.25*phim);
else
visxb = 10000;
end
 a = 0.5; b = 200; 
rhol = 2500;
Rw = 461;
T = 1353;
sigma = 0.05;
vism = 35*(10^((0.1623*y(3).^2)-0.8956*y(3)+1.2772))*(1+(10/3)*y(5));
visx = vism*visxb;
rho = rhol*(1-phib)+y(4)*phib/(Rw*T); 
u = 2500*Q/(3.14*a*b*rho);
cf(2) = (y(2)/(4*visx*u))*(y(4)-y(1)-2*sigma/y(2));
strb = (y(4)-y(1))/visx;
strb2 = (phib^(2/3))*cf(2)*(u/y(2));
phib = (y(2)^3)/(S.^3);
Machc = (y(4)/(rho*phib))^0.5;
value    = [phib > 0.5 && strb2 >1 ; (y(4)-y(1))>2400000/phib ; Machc < u ; y(1) < 100000];

isterminal = [1;1;1;1];   % Stop the integration
 
direction  = [0;0;0;0];
 
end
function cf = Tuto31(z,y,Q(k))

format long


% constants
g = 9.8; a = 0.5; b = 200;  Nb = 10^12; sigma = 0.05;
T = 1353; IFF = 10^(-13.5); kv = 4*pi/3; phim = 0.45; phimb = 0.8; 
rhol = 2500; Ro = 0.000001; B = 8.31; Rw = 461;
So = (3/(4*pi))^(1/3)*(1/Nb)^(1/3); Mw = 0.018; 
 vism = 35*(10^((0.1623*y(3).^2)-0.8956*y(3)+1.2772))*(1+(10/3)*y(5));
  tau = 20;


% variables
% dh is water diffusivity
% phib is porosity
% rho is density
% u is velocity
% n is flow index
% visx is melt-crystal viscosity
% roots are visr (vis/visx)
% vis is viscosity

dh = y(3)*exp(-8.56-19110/T);S = ((So.^3)-(Ro.^3)+(y(2).^3))^(1/3); phib = (y(2)^3)/(S.^3); rho = rhol*(1-phib)+y(4)*phib/(Rw*T); 
Q = linspace(500,510,101);
u = 2500*Q/(3.14*a*b*rho);
if y(5)<0.4409
visxb = (1-(y(5)/phim)).^(-5.25*phim);
else
visxb = 10000;
end
visx = vism*visxb;

if phib < 0.7
roots = (1-phib/phimb).^((-phimb+(5/3)*phimb*(y(2)*vism*visxb*(u/a)/sigma))/(1+(y(2)*vism*visxb*(u/a)/sigma)));
elseif phib < 0.9
roots = (1-phib/(phib+0.1)).^((-(phib+0.1)+(5/3)*(phib+0.1)*(y(2)*vism*visxb*(u/a)/sigma))/(1+(y(2)*vism*visxb*(u/a)/sigma)));
else 
roots = (1-phib).^((-1+(5/3)*(y(2)*vism*visxb*(u/a)/sigma))/(1+(y(2)*vism*visxb*(u/a)/sigma)));    
end
vis = roots*visx;
Re = 2*rho*u*a/vis;

Cr = ((-0.231+651.1/T).*((y(1)/1000000)^0.5)+(0.03424-32.57/T-0.00232)*(y(1)/1000000)); Vmo = 1/Nb; Vm = (1-y(5))*Vmo;
Machc = (y(4)/(rho*phib))^0.5;

% diff equ
cf(1) = (-rho*g-rho*(u.^2)*(24/Re+0.01)/(4*a))/(1-(u/Machc).^2);
cf(2) = (y(2)/(4*visx*u))*(y(4)-y(1)-2*sigma/y(2));
cf(6) = (1/u);
cf(5) = (1/u)*(4*kv*IFF*(y(6).^3))*exp(-kv*IFF*(y(6).^4));
cf(3) = -(4*pi*y(2)*dh)*(y(3)-Cr)/(u*Vm)+(y(3)*(cf(5)*(Vmo/Vm)));
cf(4) = (1/100)*(3*B*T*rhol*dh*(y(3)-Cr))/((y(2).^2)*Mw*u)-cf(2)*(3*y(4)/y(2));


strb2 = (phib^(2/3))*cf(2)*(u/y(2));

cf = cf(:);

Answer 1

[Star Strider]

It is difficult to follow what you are doing.

Note that in these lines in ‘Tuto31’:

dh = y(3)*exp(-8.56-19110/T);S = ((So.^3)-(Ro.^3)+(y(2).^3))^(1/3); phib = (y(2)^3)/(S.^3); rho = rhol*(1-phib)+y(4)*phib/(Rw*T); 
Q = linspace(500,510,101);              % ← Overwrites ‘Q(k)’ With This Vector
u = 2500*Q/(3.14*a*b*rho);

your over-write ‘Q’, so it never changes. (IAlso, i it supposed to be a vector or a scalar inside ‘Tuto31’?)

[Kyungmin Kim]

Thank you for your suggestion

Q is parameter (scalar) inside the "Tuto31" and I want to calculate the "Tuto31" multiple times.

changing Q start from 500, 500.1, 500.2,,,,,,,, and 510.0

and want to get the final values at each calculation.

I need to find how I can use properly for loop.

[Star Strider]

Try removing this line, in ‘Tuto31’:

Q = linspace(500,510,101);

EDIT — (12 Jan 2020 at 18:16)

I did not use your ‘events’ funciton, concentrating on the rest of the code, changing it slightly:

function cf = Tuto31(z,y,Q)
format long
% constants
g = 9.8; a = 0.5; b = 200;  Nb = 10^12; sigma = 0.05;
T = 1353; IFF = 10^(-13.5); kv = 4*pi/3; phim = 0.45; phimb = 0.8; 
rhol = 2500; Ro = 0.000001; B = 8.31; Rw = 461;
So = (3/(4*pi))^(1/3)*(1/Nb)^(1/3); Mw = 0.018; 
 vism = 35*(10^((0.1623*y(3).^2)-0.8956*y(3)+1.2772))*(1+(10/3)*y(5));
  tau = 20;
% variables
% dh is water diffusivity
% phib is porosity
% rho is density
% u is velocity
% n is flow index
% visx is melt-crystal viscosity
% roots are visr (vis/visx)
% vis is viscosity
dh = y(3)*exp(-8.56-19110/T);S = ((So.^3)-(Ro.^3)+(y(2).^3))^(1/3); phib = (y(2)^3)/(S.^3); rho = rhol*(1-phib)+y(4)*phib/(Rw*T); 
% % % Q = linspace(500,510,101);                % <— REMOVE THIS LINE
u = 2500*Q/(3.14*a*b*rho);
if y(5)<0.4409
visxb = (1-(y(5)/phim)).^(-5.25*phim);
else
visxb = 10000;
end
visx = vism*visxb;
if phib < 0.7
roots = (1-phib/phimb).^((-phimb+(5/3)*phimb*(y(2)*vism*visxb*(u/a)/sigma))/(1+(y(2)*vism*visxb*(u/a)/sigma)));
elseif phib < 0.9
roots = (1-phib/(phib+0.1)).^((-(phib+0.1)+(5/3)*(phib+0.1)*(y(2)*vism*visxb*(u/a)/sigma))/(1+(y(2)*vism*visxb*(u/a)/sigma)));
else 
roots = (1-phib).^((-1+(5/3)*(y(2)*vism*visxb*(u/a)/sigma))/(1+(y(2)*vism*visxb*(u/a)/sigma)));    
end
vis = roots*visx;
Re = 2*rho*u*a/vis;
Cr = ((-0.231+651.1/T).*((y(1)/1000000)^0.5)+(0.03424-32.57/T-0.00232)*(y(1)/1000000)); Vmo = 1/Nb; Vm = (1-y(5))*Vmo;
Machc = (y(4)/(rho*phib))^0.5;
% diff equ
cf(1) = (-rho*g-rho*(u.^2)*(24/Re+0.01)/(4*a))/(1-(u/Machc).^2);
cf(2) = (y(2)/(4*visx*u))*(y(4)-y(1)-2*sigma/y(2));
cf(6) = (1/u);
cf(5) = (1/u)*(4*kv*IFF*(y(6).^3))*exp(-kv*IFF*(y(6).^4));
cf(3) = -(4*pi*y(2)*dh)*(y(3)-Cr)/(u*Vm)+(y(3)*(cf(5)*(Vmo/Vm)));
cf(4) = (1/100)*(3*B*T*rhol*dh*(y(3)-Cr))/((y(2).^2)*Mw*u)-cf(2)*(3*y(4)/y(2));
strb2 = (phib^(2/3))*cf(2)*(u/y(2));
cf = cf(:);
end

and:

CF = cell(101,1); %assigning end value
zspan = linspace(-3200, 1000, 50); %distance start from -3200
y0 = [79000000; 0.000001; 2.85; 79100000; 0.05; 0]; % initial value 
% constants
g = 9.8; a = 0.5; b = 200; Nb = 10^12; sigma = 0.05;
T = 1353; IFF = 10^(-13.5); kv = 4*pi/3; phim = 0.6; phimb = 0.8; 
rhol = 2500; Ro = 0.000001; B = 8.31;  Rw = 461;
So = (3/(4*pi))^(1/3)*(1/Nb)^(1/3); Mw = 0.018;
% %% ode45 function
% Opt1    = odeset('Events', @myEvent1);
% Opt2    = odeset('Events', @myEvent3);
% Opt = odeset('Events', @myEvent2);
Q = linspace(500,510,11); %Changing parameter Q for each loop
for k = 1:numel(Q)
    [z,CF{k}] = ode15s(@(z,y) Tuto31(z,y,Q(k)),zspan,y0); %Assigning end value of fundction 
%     sol(k) = ode15s(@(z,y) Tuto31(z,y,Q(k)),zspan,y0,Opt); %Assigning end value of fundction 
end
% y = y(y==real(y));

figure
hold on
for k = 1:numel(Q)
    plot(z, real(CF{k}))
end
hold off
legend(compose('Q = %d\n',Q))

There is likely no need to have more than 11 elements in ‘Q’, since it does not change much, at least when plotted.

[Kyungmin Kim]

Really Thank you for your help

[Star Strider]

As always, my pleasure!