Can you help solving this differential equation using ode15i, please?

function dydz = MS_AA_imp(z, x, dx_dz ,x_p, V, Z,h)
TempK = 25 + 273.15; % Temp in Kelvin
Faraday = 96485.336;       % Faraday constant [C/mol]
R_id = 8.314;            % R constant [J/(K mol)]
VM_1 = 1.8e-5; % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS
VM_4 = 1.8e-5;          % Water molar volume [m3/mol]
D_14_inf = 6.67e-10;         % D14 value under diluted conditions and IEP [m^2/s]---substituted for D14_inf amino acid which I found on literature and used for the
%calculation of Diffusivity at different concentration(imperical relation fourn using experimental value)   D= 10^(-9) (0.665-0.393x_i)
D_24 = 1.33E-9;         % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9;         % D34 value under diluted conditions [m^2/s]
x_1 = x(1);
x_2 = x(2);
x_3 = x(3);
x_4 = x(4);
x_5 = x(5); %psi.
h=7; %hydration number at specific pH and ionic strength
%Vol_exclusion = (1 + 4*x_1 + 4*x_1^2 - 4*x_1^3 + x_1^4)/((1-x_1)^3);
hydration= x_1./((x_1).*(1-x_1.*h).*(R_id*TempK)); % new hydration formula
%CT = x_1 / VM_1 + x_2 / VM_2 + x_3 / VM_3 + x_4 / VM_4;
CT = x_1 / VM_1 + x_2 / VM_2 + x_3 / VM_3 + x_4 / VM_4;
%D_14 = D_14_inf * ( 0.21 + 0.79* exp(-4.7*x_1))/ hydration;
D_14 =10e-9 .*(0.665-0.393.*x_1)./(1+x_1.* (h./(1-x_1.* h))); % new D_14 formula
IStr = 0.5* (Z(2).^2 .* x_2/(VM_2 *1000) + Z(3).^2 .* x_3/(VM_3 *1000));
D_23 = ((D_24 + D_34)./2 .* IStr .^(0.55) ./ (abs((Z(2)) * (Z(3))).^(2.3)));
D_23V = D_23 * VM_3 *CT;
D_32V = D_23 * VM_2 *CT;
J_1 = x_p(1) * V; %think about changing to total molar flux instead of using volumetric flux V (x_tot*v=N_tot--but x_tot=1 thus N_tot=v)
J_2 = x_p(2) * V;%know that the volumetric flux here the same as the total molar flux since the concentration 
J_3 = x_p(3) * V;
J_4 = x_p(4) * V;
dpsidz = x_5;
%eq1 = Vol_exclusion * dphi_dz(1) + x_1 * Z(1) * Faraday / (R_id .* TempK)* dphi_dz(5) - (J_4 * x_1 - J_1* x_4)/(D_14 * VM_4 * CT); 
eq1 = hydration * dx_dz(1) + x_1 * Z(1) * Faraday / (R_id .* TempK)* dx_dz(5) - (J_4 * x_1 - J_1* x_4)/(D_14 * CT); 
%eq1 = hydration* diff(x_1) + x_1 * Z_AA * Faraday * diff(psi)/ ( R * TempK) == N_4 * x_1/ D_14; 
eq2 = dx_dz(2) + x_2 * Z(2) * Faraday * dx_dz(5) / ( R_id * TempK) - (J_3 * x_2 - J_2 * x_3)/ D_23 - (J_4 * x_2 - J_2 * x_4)/ (D_24 *CT);
eq3 = dx_dz(3) + x_3 * Z(3) * Faraday * dx_dz(5) / ( R_id * TempK) - (J_2 * x_3 - J_3 * x_2)/ D_23 - (J_4 * x_3 - J_3 * x_4)/ (D_34* CT);   %D_32V replaced with D_    
eq4 = x_1 + x_2 + x_3 + x_4 - 1; % contiuty
%eq5 = Z(1) * x_1/ VM_1 + Z(2)* x_2/ VM_2 + Z(3) * x_3/VM_3;
eq5 = Z(1) * x_1 + Z(2)* x_2 + Z(3) * x_3; % changed from the above line removing VM
dydz = [eq1; eq2; eq3; eq4; eq5];
end

[Torsten]

Please add your call to ode15i (tspan, initial conditions for x(1)-x(5), initialization of variables x_p, V, Z,h, call to decic etc.)

[Tadele]

function dydz = MS_AA_imp(z, x, dx_dz ,x_p, V, Z,h, phi_start)

TempK = 25 + 273.15; % Temp in Kelvin

Faraday = 96485.336; % Faraday constant [C/mol]

R_id = 8.314; % R constant [J/(K mol)]

VM_1 = 1.8e-5; % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol

VM_2 = 1.57e-5; % Na molar volume [m3/mol] * Calculated from rS

VM_3 = 4.47e-6; % Cl molar volume [m3/mol] * calcualted form rS

VM_4 = 1.8e-5; % Water molar volume [m3/mol]

D_14_inf = 6.67e-10; % D14 value under diluted conditions and IEP [m^2/s]---substituted for D14_inf amino acid which I found on literature and used for the

%calculation of Diffusivity at different concentration(imperical relation fourn using experimental value) D= 10^(-9) (0.665-0.393x_i)

D_24 = 1.33E-9; % D24 value under diluted conditions [m^2/s]

D_34 = 2.03E-9; % D34 value under diluted conditions [m^2/s]

x_1 = x(1);

x_2 = x(2);

x_3 = x(3);

x_4 = x(4);

x_5 = x(5); %psi.

V=0.000255;

h=7; %hydration number at specific pH and ionic strength

phi_start=[0 0 0];

Cp_guess_MS=[0.001 0.090 0.003];

%Vol_exclusion = (1 + 4*x_1 + 4*x_1^2 - 4*x_1^3 + x_1^4)/((1-x_1)^3);

hydration= x_1./((x_1).*(1-x_1.*h).*(R_id*TempK)); % new hydration formula

%CT = x_1 / VM_1 + x_2 / VM_2 + x_3 / VM_3 + x_4 / VM_4;

CT = x_1 / VM_1 + x_2 / VM_2 + x_3 / VM_3 + x_4 / VM_4;

%D_14 = D_14_inf * ( 0.21 + 0.79* exp(-4.7*x_1))/ hydration;

D_14 =10e-9 .*(0.665-0.393.*x_1)./(1+x_1.* (h./(1-x_1.* h))); % new D_14 formula

IStr = 0.5* (Z(2).^2 .* x_2/(VM_2 *1000) + Z(3).^2 .* x_3/(VM_3 *1000));

D_23 = ((D_24 + D_34)./2 .* IStr .^(0.55) ./ (abs((Z(2)) * (Z(3))).^(2.3)));

D_23V = D_23 * VM_3 *CT;

D_32V = D_23 * VM_2 *CT;

J_1 = x_p(1) * V; %think about changing to total molar flux instead of using volumetric flux V (x_tot*v=N_tot--but x_tot=1 thus N_tot=v)

J_2 = x_p(2) * V;%know that the volumetric flux here the same as the total molar flux since the concentration

J_3 = x_p(3) * V;

J_4 = x_p(4) * V;

dpsidz = x_5;

%eq1 = Vol_exclusion * dphi_dz(1) + x_1 * Z(1) * Faraday / (R_id .* TempK)* dphi_dz(5) - (J_4 * x_1 - J_1* x_4)/(D_14 * VM_4 * CT);

eq1 = hydration * dx_dz(1) + x_1 * Z(1) * Faraday / (R_id .* TempK)* dx_dz(5) - (J_4 * x_1 - J_1* x_4)/(D_14 * CT);

%eq1 = hydration* diff(x_1) + x_1 * Z_AA * Faraday * diff(psi)/ ( R * TempK) == N_4 * x_1/ D_14;

eq2 = dx_dz(2) + x_2 * Z(2) * Faraday * dx_dz(5) / ( R_id * TempK) - (J_3 * x_2 - J_2 * x_3)/ D_23 - (J_4 * x_2 - J_2 * x_4)/ (D_24 *CT);

eq3 = dx_dz(3) + x_3 * Z(3) * Faraday * dx_dz(5) / ( R_id * TempK) - (J_2 * x_3 - J_3 * x_2)/ D_23 - (J_4 * x_3 - J_3 * x_4)/ (D_34* CT); %D_32V replaced with D_

eq4 = x_1 + x_2 + x_3 + x_4 - 1; % contiuty

%eq5 = Z(1) * x_1/ VM_1 + Z(2)* x_2/ VM_2 + Z(3) * x_3/VM_3;

eq5 = Z(1) * x_1 + Z(2)* x_2 + Z(3) * x_3; % changed from the above line removing VM

dydz = [eq1; eq2; eq3; eq4; eq5];

[y0, yp0] = decic(@MS_AA_imp, 0, [phi_start'; 1], [1 1 0 0 0], [0 0 0 0 0], [],[],Cp_guess_MS, Jv, Z);

stepsMS = linspace(0 , z_pol, 500) ; %i think for z_pol is thickness of cp lyer

[zsol_MS, ysol_MS] = ode15i(@(z, x, dx_dz)MS_AA_imp(z, x, dx_dz ,Cp_guess_MS, Jv, Z), stepsMS, y0, yp0);

plot(zsol_MS,ysol_MS(:,1),'k:',zsol_MS,ysol_MS(:,2),zsol_MS,ysol_MS(:,3),zsol_MS,ysol_MS(:,4),'r:','LineWidth',2)

hold on

yline(0.0,'k:');

xline(delta,'r:');

xlabel("delta")

ylabel("\phi (-)")

hold off

[Tadele]

Hi Torsten,

I tried to call the ode15i based on your comment. It is on the top of this reply. Can you please help with what I missed on solving this odes?

[Torsten]

Could you arrange your code such that it can be executed ?

I.e. take the call to ode15i out of the function MS_AA_imp that is called by ode15i. Otherwise you will get an infinite loop.

[Tadele]

I get this error I don't know what the problem is.

[Torsten]

You can't run this function independently.

You must write a script in which you call ode15i and denote MS_AA_imp as the function ode15i should use to get its right-hand sides.

I suggest you start with one of the documented examples for ode15i and arrange your code accordingly.

https://de.mathworks.com/help/matlab/ref/ode15i.html

[Tadele]

function dxdz = MS_AA_imp(z,x,xp)

TempK = 25 + 273.15; % Temp in Kelvin

Faraday = 96485.336; % Faraday constant [C/mol]

R_id = 8.314; % R constant [J/(K mol)]

VM_1 = 1.8e-5; % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol

VM_2 = 1.57e-5; % Na molar volume [m3/mol] * Calculated from rS

VM_3 = 4.47e-6; % Cl molar volume [m3/mol] * calcualted form rS

VM_4 = 1.8e-5; % Water molar volume [m3/mol]

D_14_inf = 6.67e-10; % D14 value under diluted conditions and IEP [m^2/s]---substituted for D14_inf amino acid which I found on literature and used for the

%calculation of Diffusivity at different concentration(imperical relation fourn using experimental value) D= 10^(-9) (0.665-0.393x_i)

D_24 = 1.33E-9; % D24 value under diluted conditions [m^2/s]

D_34 = 2.03E-9; % D34 value under diluted conditions [m^2/s]

%x_5 = x(5); %psi.

V=0.000255;

h=7; %hydration number at specific pH and ionic strength

phi_start=[0 0 0];

Cp_guess_MS=[0.001 0.090 0.003];

%Vol_exclusion = (1 + 4*x_1 + 4*x_1^2 - 4*x_1^3 + x_1^4)/((1-x_1)^3);

hydration= x(1)/((x(1)).*(1-x(1).*h).*(R_id*TempK)); % new hydration formula

%CT = x_1 / VM_1 + x_2 / VM_2 + x_3 / VM_3 + x_4 / VM_4;

CT = x(1) / VM_1 + x(2)/ VM_2 + x(3) / VM_3 + (1-sum(x)) / VM_4;

Z=[-1 1 -1];

%D_14 = D_14_inf * ( 0.21 + 0.79* exp(-4.7*x_1))/ hydration;

D_14 =10e-9 *(0.665-0.393*x(1))/(1+x(1)* (h/(1-x(1)* h))); % new D_14 formula

IStr = 0.5* (Z(2)^2 * x(2)/(VM_2 *1000) + Z(3)^2 * x(3)/(VM_3 *1000));

D_23 = ((D_24 + D_34)/2 * IStr ^(0.55) / (abs((Z(2)) * (Z(3)))^(2.3)));

D_23V = D_23 * VM_3 *CT;

D_32V = D_23 * VM_2 *CT;

x_p=[0.001 0.002 0.0001];

x(4)=1-x(1)+x(2)+x(3);

J_1 = x_p(1) * V; %think about changing to total molar flux instead of using volumetric flux V (x_tot*v=N_tot--but x_tot=1 thus N_tot=v)

J_2 = x_p(2) * V;%know that the volumetric flux here the same as the total molar flux since the concentration

J_3 = x_p(3) * V;

J_4 = (1-sum(x_p)) * V;

%dpsidz = x(5);

dxdz = [hydration * xp(1) + x(1) * Z(1) * Faraday / (R_id .* TempK) - (J_4 * x(1) - J_1* x(4))/(D_14 * CT);

xp(2) + x(2) * Z(2) * Faraday / ( R_id * TempK) - (J_3 * x(2) - J_2 * x(3))/ D_23 - (J_4 * x(2) - J_2 * x(4))/ (D_24 *CT);

xp(3) + x(3) * Z(3) * Faraday / ( R_id * TempK) - (J_2 * x(3) - J_3 * x(2))/ D_23 - (J_4 * x(3) - J_3 * x(4))/ (D_34* CT); %D_32V replaced with D_

x(1) + x(2) + x(3) + x(4) - 1;

Z(1) * x(1) + Z(2)* x(2) + Z(3) * x(3)];

[Tadele]

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...

'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});

y0 = [1; 0; 1e-3];

yp0 = [0; 0; 0];

[y0,yp0] = decic(@MS_AA_imp,0,y0,[0 0 0],yp0,[],options);

%%

% Solve the system of DAEs using |ode15i|.

tspan = [0 4*logspace(-6,6)];

[z,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

%%

% Plot the solution components. Since the second solution component is

% small relative to the others, multiply it by |1e4| before plotting.

y(:,2) = 1e4*y(:,2);

semilogx(z,y)

ylabel('1e4 * y(:,2)')

title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

[Tadele]

Hi Torsten,

I rearrange the code in such a way that the link you have provided me shows.

The top two comments are the function and the ode calling codes. But it still shows erro when I run it. Can you please take a look at it or try it in your computer? I have four scripts that I still need to link with this ode solution so that the algorithm can run well. If this going to work, it will really help for the next step.

[Torsten]

You have 5 differential equations, but only 4 solution variables. This won't work.

What's the matter with x(5) ?

What is the z-span you want to integrate over ?

What are the initial conditions for x(1) - x(5) ?

You only copied the test example for ode15i, but didn't adapt it to your own problem.

function dxdz = MS_AA_imp(z,x,xp)
TempK = 25 + 273.15; % Temp in Kelvin
Faraday = 96485.336;       % Faraday constant [C/mol]
R_id = 8.314;            % R constant [J/(K mol)]
VM_1 = 1.8e-5; % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS
VM_4 = 1.8e-5;          % Water molar volume [m3/mol]
D_14_inf = 6.67e-10;         % D14 value under diluted conditions and IEP [m^2/s]---substituted for D14_inf amino acid which I found on literature and used for the
%calculation of Diffusivity at different concentration(imperical relation fourn using experimental value)   D= 10^(-9) (0.665-0.393x_i)
D_24 = 1.33E-9;         % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9;         % D34 value under diluted conditions [m^2/s]
%x_5 = x(5); %psi.
V=0.000255;
h=7; %hydration number at specific pH and ionic strength
phi_start=[0 0 0];
Cp_guess_MS=[0.001 0.090 0.003];
%Vol_exclusion = (1 + 4*x_1 + 4*x_1^2 - 4*x_1^3 + x_1^4)/((1-x_1)^3);
hydration= x(1)/((x(1)).*(1-x(1).*h).*(R_id*TempK)); % new hydration formula
%CT = x_1 / VM_1 + x_2 / VM_2 + x_3 / VM_3 + x_4 / VM_4;
CT = x(1) / VM_1 + x(2)/ VM_2 + x(3) / VM_3 + (1-sum(x)) / VM_4;
Z=[-1 1 -1];
%D_14 = D_14_inf * ( 0.21 + 0.79* exp(-4.7*x_1))/ hydration;
D_14 =10e-9 *(0.665-0.393*x(1))/(1+x(1)* (h/(1-x(1)* h))); % new D_14 formula
IStr = 0.5* (Z(2)^2 * x(2)/(VM_2 *1000) + Z(3)^2 * x(3)/(VM_3 *1000));
D_23 = ((D_24 + D_34)/2 * IStr ^(0.55) / (abs((Z(2)) * (Z(3)))^(2.3)));
D_23V = D_23 * VM_3 *CT;
D_32V = D_23 * VM_2 *CT;
x_p=[0.001 0.002 0.0001];
x(4)=1-x(1)+x(2)+x(3);
J_1 = x_p(1) * V; %think about changing to total molar flux instead of using volumetric flux V (x_tot*v=N_tot--but x_tot=1 thus N_tot=v)
J_2 = x_p(2) * V;%know that the volumetric flux here the same as the total molar flux since the concentration 
J_3 = x_p(3) * V;
J_4 = (1-sum(x_p)) * V;
%dpsidz = x(5);
dxdz = [hydration * xp(1) + x(1) * Z(1) * Faraday / (R_id .* TempK) - (J_4 * x(1) - J_1* x(4))/(D_14 * CT); 
 xp(2) + x(2) * Z(2) * Faraday  / ( R_id * TempK) - (J_3 * x(2) - J_2 * x(3))/ D_23 - (J_4 * x(2) - J_2 * x(4))/ (D_24 *CT);
 xp(3) + x(3) * Z(3) * Faraday  / ( R_id * TempK) - (J_2 * x(3) - J_3 * x(2))/ D_23 - (J_4 * x(3) - J_3 * x(4))/ (D_34* CT);   %D_32V replaced with D_    
 x(1) + x(2) + x(3) + x(4) - 1; 
 Z(1) * x(1) + Z(2)* x(2) + Z(3) * x(3)];
end

[Tadele]

I answered some of the questions you raised by editing the text of your comment. Below are the answers.

You have 5 differential equations, but only 4 solution variables. This won't work.

What's the matter with x(5) ? --->>this is anther differential equation with different variable. But I wanted to calculate it as a fitting paramets that will be calculated automatically when othe variables(x1-x4) are solved.

What is the z-span you want to integrate over ?-------->> Z span is from zero to 1.23e-5

What are the initial conditions for x(1) - x(5) ? -------->>the initrial conditions are lab concentrations that will be varied for different experiments. But for now it can be used as x=[ 0.001 0.002 0.003 0.994 0]

Can you help where I should use them?

[Torsten]

What's the matter with x(5) ? --->>this is anther differential equation with different variable. But I wanted to calculate it as a fitting paramets that will be calculated automatically when othe variables(x1-x4) are solved.

How ?

In the testcode above:

y0 are your initial conditions.

tspan is Z span.

Adapt the call to decic to your problem now.

[Tadele]

about x(5), I removed it form the first script I posted to make it easier for someone who helps solving the equations, hoping that I will add that after I get help solving the equations without varieble x(5). But it becomes difficult now to solve it for me.

I also tried the test code you have sent me by adding the intial conditions and tsapan to adapt them to the decic. It is still showing me error. Does it work for you (if you have tried it)? Otherwise is there any other way that you recommend me to solve it?

If you have time can you also suggest me how to use the output of one function as an input for the other function so that I can run an algorithm that converges a guess value and calculated value? If it is easier with codes, I can post it here.

Answer 1

[Jan]

dx_dz = ???
x_p = ???
V =  ???
Z = ???
h = ???
tspan = [0, 10]; % ???
x0 = rand(1, 5); % ???
[Z, Y] = ode15s(@(z, x) MS_AA_imp(z, x, dx_dz, x_p, V, Z, h), tspan, x0);

Although the names of the variables do not matter, "dx_dz", "dydz" and "z", "x" is confusing.

[Torsten]

@Jan

dx_dz are input to "MS_AA_imp" from ODE15I and are the time-derivatives of the variables solved for.

Since the ODE system does not seem to be explicit in dx_dz, it will at least be necessary to use ODE15S together with the mass matrix option.

[Jan]

@Torsten: I know. I've mentioned the names of the variables due to:

function dydz = MS_AA_imp(z, x, ...
%         ^                  ^

This is not an error, but maybe confusing, especially if "dx_dz" is occurring also, even if it has a different meaning.