Good Evening, Im a physcis student, but now i want to study on geoscience and particularly using MATLAB to help me to integrating all my physics background into analyzing, visualizing, etc. with geoscience case. Is there any advice where to start into this field? thx in advanceGood Evening, Im a physcis student, but now i want to study on geoscience and particularly using MATLAB to help me to integrating all my physics background into analyzing, visualizing, etc. with geoscience case. Is there any advice where to start into this field? thx in advance Good Evening, Im a physcis student, but now i want to study on geoscience and particularly using MATLAB to help me to integrating all my physics background into analyzing, visualizing, etc. with geoscience case. Is there any advice where to start into this field? thx in advance geoscience, physics, mapping toolbox, matlab MATLAB Answers — New Questions

​

I need to project 3D nodes (nodes) onto a 2D plane that I got from the code in FEX.
load xyz

fobjPlane=planarFit(xyz);
for i=1:2
subplot(1,2,i);
[hL,hD]=plot(fobjPlane); %Visualize the fit
end
legend([hL,hD],’Plane fit’, ‘XYZ samples’, ‘Location’,’northoutside’,’FontSize’,15);
subplot(1,2,1)
view([-69.5 -24.6])
subplot(1,2,2)
view([9.6 -35.2])
So I need to always get nodes in space however that lie on that plan.I need to project 3D nodes (nodes) onto a 2D plane that I got from the code in FEX.
load xyz

fobjPlane=planarFit(xyz);
for i=1:2
subplot(1,2,i);
[hL,hD]=plot(fobjPlane); %Visualize the fit
end
legend([hL,hD],’Plane fit’, ‘XYZ samples’, ‘Location’,’northoutside’,’FontSize’,15);
subplot(1,2,1)
view([-69.5 -24.6])
subplot(1,2,2)
view([9.6 -35.2])
So I need to always get nodes in space however that lie on that plan. I need to project 3D nodes (nodes) onto a 2D plane that I got from the code in FEX.
load xyz

fobjPlane=planarFit(xyz);
for i=1:2
subplot(1,2,i);
[hL,hD]=plot(fobjPlane); %Visualize the fit
end
legend([hL,hD],’Plane fit’, ‘XYZ samples’, ‘Location’,’northoutside’,’FontSize’,15);
subplot(1,2,1)
view([-69.5 -24.6])
subplot(1,2,2)
view([9.6 -35.2])
So I need to always get nodes in space however that lie on that plan. nodes, projection, orthogonal, plan MATLAB Answers — New Questions

​

I have make a code that do the scatter points in a graph, but i need to find the number of scatter points of each color? i have attached the sample code and graph for explaintion
for h=num_fixed_nodes+1:1:num_nodes
%Generate random nodes within the region
S1(h).xd = position_region(i,1) + rand(1,1)*region_width;
S1(h).yd = position_region(i,2) + rand(1,1)*region_height;
S1(h).G=0;
S1(h).id=h;
S1(h).type=’N’;
S1(h).temp = interp2(aa1,bb1,temp_values_1,S1(h).xd,S1(h).yd);
S1(h).E=Eo*rand(1,1);
Et=Et+S1(h).E;
S1(h).node_status = S1(h).temp<thresh_temp;
if(S1(h).node_status==1)
scatter(S1(h).xd,S1(h).yd, ‘filled’,’MarkerFaceColor’,’g’);
else
scatter(S1(h).xd,S1(h).yd, ‘filled’,’MarkerFaceColor’,’r’);
end
endI have make a code that do the scatter points in a graph, but i need to find the number of scatter points of each color? i have attached the sample code and graph for explaintion
for h=num_fixed_nodes+1:1:num_nodes
%Generate random nodes within the region
S1(h).xd = position_region(i,1) + rand(1,1)*region_width;
S1(h).yd = position_region(i,2) + rand(1,1)*region_height;
S1(h).G=0;
S1(h).id=h;
S1(h).type=’N’;
S1(h).temp = interp2(aa1,bb1,temp_values_1,S1(h).xd,S1(h).yd);
S1(h).E=Eo*rand(1,1);
Et=Et+S1(h).E;
S1(h).node_status = S1(h).temp<thresh_temp;
if(S1(h).node_status==1)
scatter(S1(h).xd,S1(h).yd, ‘filled’,’MarkerFaceColor’,’g’);
else
scatter(S1(h).xd,S1(h).yd, ‘filled’,’MarkerFaceColor’,’r’);
end
end I have make a code that do the scatter points in a graph, but i need to find the number of scatter points of each color? i have attached the sample code and graph for explaintion
for h=num_fixed_nodes+1:1:num_nodes
%Generate random nodes within the region
S1(h).xd = position_region(i,1) + rand(1,1)*region_width;
S1(h).yd = position_region(i,2) + rand(1,1)*region_height;
S1(h).G=0;
S1(h).id=h;
S1(h).type=’N’;
S1(h).temp = interp2(aa1,bb1,temp_values_1,S1(h).xd,S1(h).yd);
S1(h).E=Eo*rand(1,1);
Et=Et+S1(h).E;
S1(h).node_status = S1(h).temp<thresh_temp;
if(S1(h).node_status==1)
scatter(S1(h).xd,S1(h).yd, ‘filled’,’MarkerFaceColor’,’g’);
else
scatter(S1(h).xd,S1(h).yd, ‘filled’,’MarkerFaceColor’,’r’);
end
end plotting, sum MATLAB Answers — New Questions

​

I want to use some python library to create a PGM for joint-distribution modeling and probabilistic inference. To be more specific, I have a time series data, which consists of many time slices. There are multiple variables in each time slices. Some of the variables are discrete, others are continuous. I want to do structure learning and parameter learning with this data and finally create a PGM. But I couldn’t find any suitable python library for that.
I found a library named PyPGM,

But it seems that it doesn’t support time series and hybrid(discrete and continuous)I want to use some python library to create a PGM for joint-distribution modeling and probabilistic inference. To be more specific, I have a time series data, which consists of many time slices. There are multiple variables in each time slices. Some of the variables are discrete, others are continuous. I want to do structure learning and parameter learning with this data and finally create a PGM. But I couldn’t find any suitable python library for that.
I found a library named PyPGM,

But it seems that it doesn’t support time series and hybrid(discrete and continuous) I want to use some python library to create a PGM for joint-distribution modeling and probabilistic inference. To be more specific, I have a time series data, which consists of many time slices. There are multiple variables in each time slices. Some of the variables are discrete, others are continuous. I want to do structure learning and parameter learning with this data and finally create a PGM. But I couldn’t find any suitable python library for that.
I found a library named PyPGM,

But it seems that it doesn’t support time series and hybrid(discrete and continuous) python MATLAB Answers — New Questions

​

Hello everyone!! I’m trying to plot: x,y,X,Y. But when my matlab try to plot y, the following error appears: Error using plot Vectors must be the same length. Error in doc (lne 94) plot(t, y);

Please fell free to check the code below:
Fs = 5000;
T = 1/Fs;
t = 0:T:1-T;
NFFT = 2^nextpow2(length(t));
x = cos(2*pi*100*t) + cos(2*pi*500*t) + cos(2*pi*1000*t);
[X, f, ~] = fftm(x, Fs, NFFT);
H = double(f <= 600);
Y = X .* H;
y = ifft(Y);
figure;
subplot(2,1,1);
plot(t, x);
title(‘x(t)’);
xlabel(‘Time(s)’);
ylabel(‘Amp’);
subplot(2,1,2);
plot(t, y);
title(‘y(t)’);
xlabel(‘Time(s)’);
ylabel(‘Amp’);
figure;
subplot(2,1,1);
plot(f, abs(X));
title(‘X(omega)’);
xlabel(‘Freq(Hz)’);
ylabel(‘Mag’);
subplot(2,1,2);
plot(f, abs(Y));
title(‘Y(omega)’);
xlabel(‘Freq(Hz)’);
ylabel(‘Mag’);

Oh, the following function fftm in line 6, is a function downloaded by me, it’s a fft function modified, and you can see it below:

% Modified fft algorithm
%
% This function estimates the fft of a signal x using the Fs sampling
% frequency, NFFT bins, and returns the frequencies values along with the spectrum X.
%
% [X,f] = fftm(x,Fs,NFFT)
%
% Example:
% Fs = 1e3;
% t = 0:0.001:1-0.001;
% x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
% [X,f]=fftm(x,Fs,2000);
% plot(f,abs(X));
% http://www.mathworks.com/help/signal/ug/amplitude-estimation-and-zero-padding.html
% http://www.mathworks.com/help/signal/ug/psd-estimate-using-fft.html
function [xdft,f, psdx] = fftm(x,Fs,NFFT)
L = length(x);
if nargin < 3
NFFT = 2^nextpow2(L);
end
xdft = fft(x,NFFT);
xdft = xdft(1:length(xdft)/2+1);
psdx = (1/(Fs*L)).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
psdx = 10*log10(psdx);
xdft = 1/L.*xdft;
xdft(2:end-1) = 2*xdft(2:end-1);
f = 0:Fs/NFFT:Fs/2;Hello everyone!! I’m trying to plot: x,y,X,Y. But when my matlab try to plot y, the following error appears: Error using plot Vectors must be the same length. Error in doc (lne 94) plot(t, y);

Please fell free to check the code below:
Fs = 5000;
T = 1/Fs;
t = 0:T:1-T;
NFFT = 2^nextpow2(length(t));
x = cos(2*pi*100*t) + cos(2*pi*500*t) + cos(2*pi*1000*t);
[X, f, ~] = fftm(x, Fs, NFFT);
H = double(f <= 600);
Y = X .* H;
y = ifft(Y);
figure;
subplot(2,1,1);
plot(t, x);
title(‘x(t)’);
xlabel(‘Time(s)’);
ylabel(‘Amp’);
subplot(2,1,2);
plot(t, y);
title(‘y(t)’);
xlabel(‘Time(s)’);
ylabel(‘Amp’);
figure;
subplot(2,1,1);
plot(f, abs(X));
title(‘X(omega)’);
xlabel(‘Freq(Hz)’);
ylabel(‘Mag’);
subplot(2,1,2);
plot(f, abs(Y));
title(‘Y(omega)’);
xlabel(‘Freq(Hz)’);
ylabel(‘Mag’);

Oh, the following function fftm in line 6, is a function downloaded by me, it’s a fft function modified, and you can see it below:

% Modified fft algorithm
%
% This function estimates the fft of a signal x using the Fs sampling
% frequency, NFFT bins, and returns the frequencies values along with the spectrum X.
%
% [X,f] = fftm(x,Fs,NFFT)
%
% Example:
% Fs = 1e3;
% t = 0:0.001:1-0.001;
% x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
% [X,f]=fftm(x,Fs,2000);
% plot(f,abs(X));
% http://www.mathworks.com/help/signal/ug/amplitude-estimation-and-zero-padding.html
% http://www.mathworks.com/help/signal/ug/psd-estimate-using-fft.html
function [xdft,f, psdx] = fftm(x,Fs,NFFT)
L = length(x);
if nargin < 3
NFFT = 2^nextpow2(L);
end
xdft = fft(x,NFFT);
xdft = xdft(1:length(xdft)/2+1);
psdx = (1/(Fs*L)).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
psdx = 10*log10(psdx);
xdft = 1/L.*xdft;
xdft(2:end-1) = 2*xdft(2:end-1);
f = 0:Fs/NFFT:Fs/2; Hello everyone!! I’m trying to plot: x,y,X,Y. But when my matlab try to plot y, the following error appears: Error using plot Vectors must be the same length. Error in doc (lne 94) plot(t, y);

Please fell free to check the code below:
Fs = 5000;
T = 1/Fs;
t = 0:T:1-T;
NFFT = 2^nextpow2(length(t));
x = cos(2*pi*100*t) + cos(2*pi*500*t) + cos(2*pi*1000*t);
[X, f, ~] = fftm(x, Fs, NFFT);
H = double(f <= 600);
Y = X .* H;
y = ifft(Y);
figure;
subplot(2,1,1);
plot(t, x);
title(‘x(t)’);
xlabel(‘Time(s)’);
ylabel(‘Amp’);
subplot(2,1,2);
plot(t, y);
title(‘y(t)’);
xlabel(‘Time(s)’);
ylabel(‘Amp’);
figure;
subplot(2,1,1);
plot(f, abs(X));
title(‘X(omega)’);
xlabel(‘Freq(Hz)’);
ylabel(‘Mag’);
subplot(2,1,2);
plot(f, abs(Y));
title(‘Y(omega)’);
xlabel(‘Freq(Hz)’);
ylabel(‘Mag’);

Oh, the following function fftm in line 6, is a function downloaded by me, it’s a fft function modified, and you can see it below:

% Modified fft algorithm
%
% This function estimates the fft of a signal x using the Fs sampling
% frequency, NFFT bins, and returns the frequencies values along with the spectrum X.
%
% [X,f] = fftm(x,Fs,NFFT)
%
% Example:
% Fs = 1e3;
% t = 0:0.001:1-0.001;
% x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
% [X,f]=fftm(x,Fs,2000);
% plot(f,abs(X));
% http://www.mathworks.com/help/signal/ug/amplitude-estimation-and-zero-padding.html
% http://www.mathworks.com/help/signal/ug/psd-estimate-using-fft.html
function [xdft,f, psdx] = fftm(x,Fs,NFFT)
L = length(x);
if nargin < 3
NFFT = 2^nextpow2(L);
end
xdft = fft(x,NFFT);
xdft = xdft(1:length(xdft)/2+1);
psdx = (1/(Fs*L)).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
psdx = 10*log10(psdx);
xdft = 1/L.*xdft;
xdft(2:end-1) = 2*xdft(2:end-1);
f = 0:Fs/NFFT:Fs/2; vector, vectors must be the same length, error, matlab, vectors error, plot, plot error, fourier MATLAB Answers — New Questions

​

Hi!
I’m trying to model gas separation of a gas with 4 components (H2, N2, N3H and O2) using a hollow fiber membrane. In papers I found the following equations:

(the change in the retentate flow rate)
( the change in the mole fraction of component i in the retentate)
( the change in the mole fraction of component i in the permeate)
I am solving these equations by first using the backwards finite difference method which turns the first two equations into:

Next, I am trying to use Newton’s method to solve all 9 equations together. This is my code:
%% Gas separation in a cocurrent hollow fiber membrane
thickness_membrane = 1e-6; % thickness membrane [m]
Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_O2 = 37*3.35e-16;
R=8.314;
T = 273.15+25; % correlation temperature [K]
Per_H2 = Perm_H2/thickness_membrane; % Permeance of H2 [mol/m2 s Pa]
Per_N2 = Perm_N2/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_NH3 = Perm_NH3/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_O2 = Perm_O2/thickness_membrane;

%% Input parameters
F_feed = 1133.95/3.6; % feed [mol/s]
x_H2_F = 0.0268; % [-]
x_N2_F =0.9682; % [-]
x_NH3_F = 0.0049; % [-]
x_O2_F = 0.0001;

P_F = 35e5; % pressure feed [Pa]
P_P = 1e5; % pressure Permeate [Pa]

%% Assumptions membrane
L_fiber = 10;
r_fiber = (0.75e-3)/2;
n_fiber = 4000;

% Define the mesh
step= 80;
mesh = 0:(L_fiber/step):L_fiber; % linear space vector
nmesh= numel(mesh);
h = mesh(2)-mesh(1);
visc_H2 = 88e-6;
visc_N2 = 17.82e-6; % [Pa s]
visc_NH3 = 9.9e-6;
visc_O2 = 20.4e-6;

mu = visc_N2; % most present

syms xH2 yH2 Ff xN2 yN2 xNH3 yNH3 xO2 yO2
variables = [ Ff; xH2; xN2 ; xNH3; xO2; yH2; yN2; yNH3; yO2];
Final_results = zeros(10,10,numel(mesh));
xH2prev = x_H2_F;
xN2prev = x_N2_F;
xNH3prev = x_NH3_F;
xO2prev = x_O2_F;
%P_Pprev= P_P;
Ffprev = F_feed;

BETA = P_P/P_F;
for z = 1:nmesh
%% Equations to solve making use of the backward difference method
eq1 = (Ff-Ffprev)/h + 2*3.14*r_fiber*L_fiber*n_fiber*((Per_H2*(xH2*P_F + yH2*P_P)) + (Per_NH3*(xNH3*P_F + yNH3*P_P)) + (Per_N2*(xN2*P_F + yN2*P_P)) + (Per_O2*(xO2*P_F + yO2*P_P)) );

eq2 = (xH2- xH2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_H2*(xH2*P_F + yH2*P_P) + xH2*(Ff-Ffprev)/h);
eq3 = (xN2- xN2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_N2*(xN2*P_F + yN2*P_P) + xN2*(Ff-Ffprev)/h);
eq4 = (xNH3- xNH3prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_NH3*(xNH3*P_F + yNH3*P_P) + xNH3*(Ff-Ffprev)/h);
eq5 = (xO2- xO2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_O2*(xO2*P_F + yO2*P_P) + xO2*(Ff-Ffprev)/h);

eq6 = yH2 – (Per_H2*xH2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_H2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq7 = yN2 – (Per_N2*xN2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_N2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq8 = yNH3 – (Per_NH3*xNH3 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_NH3*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq9 = yO2 – (Per_O2*xO2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_O2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));

%eq10 = ((Pp – Ppprev)/h) -( (8*R*T*mu*(F_feed-Ff))/ (3.14*r_fiber^4*n_fiber*Pp));

F = [eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9];
J=jacobian([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9],variables);
x_0 = [10 ;0.5; 0.001; 0.1; 0.0001; 0.3; 0.01; 1; 0.001];

Final_results(1,1:9,1) = x_0′;
iterations = 100;
for iter=1:iterations
%% Newton
% Substitute the initial values into the Jacobian matrix
J_subs = subs(J, variables, x_0);
F_subs = subs(F,variables,x_0);
aug_matrix = [J_subs, -F_subs];
n= 9;
for i = 1:n-1
for j = i+1:n
factor = aug_matrix(j, i) / aug_matrix(i, i);
aug_matrix(j, 🙂 = aug_matrix(j, 🙂 – factor * aug_matrix(i, :);
end
end

y_0 = zeros(n, 1);
y_0(n) = aug_matrix(n, n+1) / aug_matrix(n,n);
for i = n-1:-1:1
y_0(i) = (aug_matrix(i, n+1) – aug_matrix(i, i+1:n) * y_0(i+1:n)) / aug_matrix(i, i);
end

x_result = y_0 + x_0;
Final_results(iter+1,1:9,z) = x_result’;
err = norm(x_result – x_0);
Final_results(iter+1,10,z) = err;
x_0 = x_result;
if err <= 1e-8

%P_Pprev = x_result(10);

disp([‘Converged after ‘, num2str(iter), ‘ iterations’]);
break;
end
end
xH2prev = x_0(2);
xN2prev = x_0(3);
xNH3prev = x_0(4);
xO2prev = x_0(5);
Ffprev = x_0(1);

end

solution_Ff_Per_step = [];
solution_xH2_Per_step = [];
solution_xN2_Per_step = [];
solution_xNH3_Per_step = [];
solution_xO2_Per_step = [];
solution_yH2_Per_step = [];
solution_yN2_Per_step = [];
solution_yNH3_Per_step = [];
solution_yO2_Per_step = [];
solution_P_P_Per_step = [];

for m=1:numel(mesh)
solution_Ff_Per_step = [Final_results(end,1,m) solution_Ff_Per_step];
solution_xH2_Per_step = [Final_results(end,2,m) solution_xH2_Per_step];
solution_xN2_Per_step = [Final_results(end,3,m) solution_xN2_Per_step];
solution_xNH3_Per_step = [Final_results(end,4,m) solution_xNH3_Per_step];
solution_xO2_Per_step = [Final_results(end,5,m) solution_xO2_Per_step];
solution_yH2_Per_step = [Final_results(end,6,m) solution_yH2_Per_step];
solution_yN2_Per_step = [Final_results(end,7,m) solution_yN2_Per_step];
solution_yNH3_Per_step = [Final_results(end,8,m) solution_yNH3_Per_step];
solution_yO2_Per_step = [Final_results(end,9,m) solution_yO2_Per_step];
%solution_P_P_Per_step = [Final_results(end,10,m) solution_P_P_Per_step];

end
figure(1)
plot(mesh,solution_Ff_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Flow rate retentate [m/s]’)
figure(2)
plot(mesh,solution_xH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in retentate’)
figure(3)
plot(mesh,solution_xN2_Per_step )
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in retentate’)
figure(4)
plot(mesh,solution_xNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in retentate’)
figure(5)
plot(mesh,solution_xO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in retentate’)
figure(6)
plot(mesh,solution_yH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in Permeate’)
figure(7)
plot(mesh,solution_yN2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in Permeate’)
figure(8)
plot(mesh,solution_yNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in Permeate’)
figure(9)
plot(mesh,solution_yO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in Permeate’)
figure(10)
% plot(mesh,solution_P_P_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Permeate pressure [bar]’)
For some reason, when the method converges it will always return the initial values of xH2prev, xNH3prev, xO2prev. xN2prev and Ffprev, meaning that there is no change in mole fractions along the membrane. When I try to change the initial conditions, the graphs show spikes in random positions whilst being 0 in other positions along the membrane. This should not be the case. Could you please help me figure out what went wrong in my code?Hi!
I’m trying to model gas separation of a gas with 4 components (H2, N2, N3H and O2) using a hollow fiber membrane. In papers I found the following equations:

(the change in the retentate flow rate)
( the change in the mole fraction of component i in the retentate)
( the change in the mole fraction of component i in the permeate)
I am solving these equations by first using the backwards finite difference method which turns the first two equations into:

Next, I am trying to use Newton’s method to solve all 9 equations together. This is my code:
%% Gas separation in a cocurrent hollow fiber membrane
thickness_membrane = 1e-6; % thickness membrane [m]
Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_O2 = 37*3.35e-16;
R=8.314;
T = 273.15+25; % correlation temperature [K]
Per_H2 = Perm_H2/thickness_membrane; % Permeance of H2 [mol/m2 s Pa]
Per_N2 = Perm_N2/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_NH3 = Perm_NH3/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_O2 = Perm_O2/thickness_membrane;

%% Input parameters
F_feed = 1133.95/3.6; % feed [mol/s]
x_H2_F = 0.0268; % [-]
x_N2_F =0.9682; % [-]
x_NH3_F = 0.0049; % [-]
x_O2_F = 0.0001;

P_F = 35e5; % pressure feed [Pa]
P_P = 1e5; % pressure Permeate [Pa]

%% Assumptions membrane
L_fiber = 10;
r_fiber = (0.75e-3)/2;
n_fiber = 4000;

% Define the mesh
step= 80;
mesh = 0:(L_fiber/step):L_fiber; % linear space vector
nmesh= numel(mesh);
h = mesh(2)-mesh(1);
visc_H2 = 88e-6;
visc_N2 = 17.82e-6; % [Pa s]
visc_NH3 = 9.9e-6;
visc_O2 = 20.4e-6;

mu = visc_N2; % most present

syms xH2 yH2 Ff xN2 yN2 xNH3 yNH3 xO2 yO2
variables = [ Ff; xH2; xN2 ; xNH3; xO2; yH2; yN2; yNH3; yO2];
Final_results = zeros(10,10,numel(mesh));
xH2prev = x_H2_F;
xN2prev = x_N2_F;
xNH3prev = x_NH3_F;
xO2prev = x_O2_F;
%P_Pprev= P_P;
Ffprev = F_feed;

BETA = P_P/P_F;
for z = 1:nmesh
%% Equations to solve making use of the backward difference method
eq1 = (Ff-Ffprev)/h + 2*3.14*r_fiber*L_fiber*n_fiber*((Per_H2*(xH2*P_F + yH2*P_P)) + (Per_NH3*(xNH3*P_F + yNH3*P_P)) + (Per_N2*(xN2*P_F + yN2*P_P)) + (Per_O2*(xO2*P_F + yO2*P_P)) );

eq2 = (xH2- xH2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_H2*(xH2*P_F + yH2*P_P) + xH2*(Ff-Ffprev)/h);
eq3 = (xN2- xN2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_N2*(xN2*P_F + yN2*P_P) + xN2*(Ff-Ffprev)/h);
eq4 = (xNH3- xNH3prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_NH3*(xNH3*P_F + yNH3*P_P) + xNH3*(Ff-Ffprev)/h);
eq5 = (xO2- xO2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_O2*(xO2*P_F + yO2*P_P) + xO2*(Ff-Ffprev)/h);

eq6 = yH2 – (Per_H2*xH2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_H2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq7 = yN2 – (Per_N2*xN2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_N2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq8 = yNH3 – (Per_NH3*xNH3 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_NH3*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq9 = yO2 – (Per_O2*xO2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_O2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));

%eq10 = ((Pp – Ppprev)/h) -( (8*R*T*mu*(F_feed-Ff))/ (3.14*r_fiber^4*n_fiber*Pp));

F = [eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9];
J=jacobian([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9],variables);
x_0 = [10 ;0.5; 0.001; 0.1; 0.0001; 0.3; 0.01; 1; 0.001];

Final_results(1,1:9,1) = x_0′;
iterations = 100;
for iter=1:iterations
%% Newton
% Substitute the initial values into the Jacobian matrix
J_subs = subs(J, variables, x_0);
F_subs = subs(F,variables,x_0);
aug_matrix = [J_subs, -F_subs];
n= 9;
for i = 1:n-1
for j = i+1:n
factor = aug_matrix(j, i) / aug_matrix(i, i);
aug_matrix(j, 🙂 = aug_matrix(j, 🙂 – factor * aug_matrix(i, :);
end
end

y_0 = zeros(n, 1);
y_0(n) = aug_matrix(n, n+1) / aug_matrix(n,n);
for i = n-1:-1:1
y_0(i) = (aug_matrix(i, n+1) – aug_matrix(i, i+1:n) * y_0(i+1:n)) / aug_matrix(i, i);
end

x_result = y_0 + x_0;
Final_results(iter+1,1:9,z) = x_result’;
err = norm(x_result – x_0);
Final_results(iter+1,10,z) = err;
x_0 = x_result;
if err <= 1e-8

%P_Pprev = x_result(10);

disp([‘Converged after ‘, num2str(iter), ‘ iterations’]);
break;
end
end
xH2prev = x_0(2);
xN2prev = x_0(3);
xNH3prev = x_0(4);
xO2prev = x_0(5);
Ffprev = x_0(1);

end

solution_Ff_Per_step = [];
solution_xH2_Per_step = [];
solution_xN2_Per_step = [];
solution_xNH3_Per_step = [];
solution_xO2_Per_step = [];
solution_yH2_Per_step = [];
solution_yN2_Per_step = [];
solution_yNH3_Per_step = [];
solution_yO2_Per_step = [];
solution_P_P_Per_step = [];

for m=1:numel(mesh)
solution_Ff_Per_step = [Final_results(end,1,m) solution_Ff_Per_step];
solution_xH2_Per_step = [Final_results(end,2,m) solution_xH2_Per_step];
solution_xN2_Per_step = [Final_results(end,3,m) solution_xN2_Per_step];
solution_xNH3_Per_step = [Final_results(end,4,m) solution_xNH3_Per_step];
solution_xO2_Per_step = [Final_results(end,5,m) solution_xO2_Per_step];
solution_yH2_Per_step = [Final_results(end,6,m) solution_yH2_Per_step];
solution_yN2_Per_step = [Final_results(end,7,m) solution_yN2_Per_step];
solution_yNH3_Per_step = [Final_results(end,8,m) solution_yNH3_Per_step];
solution_yO2_Per_step = [Final_results(end,9,m) solution_yO2_Per_step];
%solution_P_P_Per_step = [Final_results(end,10,m) solution_P_P_Per_step];

end
figure(1)
plot(mesh,solution_Ff_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Flow rate retentate [m/s]’)
figure(2)
plot(mesh,solution_xH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in retentate’)
figure(3)
plot(mesh,solution_xN2_Per_step )
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in retentate’)
figure(4)
plot(mesh,solution_xNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in retentate’)
figure(5)
plot(mesh,solution_xO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in retentate’)
figure(6)
plot(mesh,solution_yH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in Permeate’)
figure(7)
plot(mesh,solution_yN2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in Permeate’)
figure(8)
plot(mesh,solution_yNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in Permeate’)
figure(9)
plot(mesh,solution_yO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in Permeate’)
figure(10)
% plot(mesh,solution_P_P_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Permeate pressure [bar]’)
For some reason, when the method converges it will always return the initial values of xH2prev, xNH3prev, xO2prev. xN2prev and Ffprev, meaning that there is no change in mole fractions along the membrane. When I try to change the initial conditions, the graphs show spikes in random positions whilst being 0 in other positions along the membrane. This should not be the case. Could you please help me figure out what went wrong in my code? Hi!
I’m trying to model gas separation of a gas with 4 components (H2, N2, N3H and O2) using a hollow fiber membrane. In papers I found the following equations:

(the change in the retentate flow rate)
( the change in the mole fraction of component i in the retentate)
( the change in the mole fraction of component i in the permeate)
I am solving these equations by first using the backwards finite difference method which turns the first two equations into:

Next, I am trying to use Newton’s method to solve all 9 equations together. This is my code:
%% Gas separation in a cocurrent hollow fiber membrane
thickness_membrane = 1e-6; % thickness membrane [m]
Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_O2 = 37*3.35e-16;
R=8.314;
T = 273.15+25; % correlation temperature [K]
Per_H2 = Perm_H2/thickness_membrane; % Permeance of H2 [mol/m2 s Pa]
Per_N2 = Perm_N2/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_NH3 = Perm_NH3/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_O2 = Perm_O2/thickness_membrane;

%% Input parameters
F_feed = 1133.95/3.6; % feed [mol/s]
x_H2_F = 0.0268; % [-]
x_N2_F =0.9682; % [-]
x_NH3_F = 0.0049; % [-]
x_O2_F = 0.0001;

P_F = 35e5; % pressure feed [Pa]
P_P = 1e5; % pressure Permeate [Pa]

%% Assumptions membrane
L_fiber = 10;
r_fiber = (0.75e-3)/2;
n_fiber = 4000;

% Define the mesh
step= 80;
mesh = 0:(L_fiber/step):L_fiber; % linear space vector
nmesh= numel(mesh);
h = mesh(2)-mesh(1);
visc_H2 = 88e-6;
visc_N2 = 17.82e-6; % [Pa s]
visc_NH3 = 9.9e-6;
visc_O2 = 20.4e-6;

mu = visc_N2; % most present

syms xH2 yH2 Ff xN2 yN2 xNH3 yNH3 xO2 yO2
variables = [ Ff; xH2; xN2 ; xNH3; xO2; yH2; yN2; yNH3; yO2];
Final_results = zeros(10,10,numel(mesh));
xH2prev = x_H2_F;
xN2prev = x_N2_F;
xNH3prev = x_NH3_F;
xO2prev = x_O2_F;
%P_Pprev= P_P;
Ffprev = F_feed;

BETA = P_P/P_F;
for z = 1:nmesh
%% Equations to solve making use of the backward difference method
eq1 = (Ff-Ffprev)/h + 2*3.14*r_fiber*L_fiber*n_fiber*((Per_H2*(xH2*P_F + yH2*P_P)) + (Per_NH3*(xNH3*P_F + yNH3*P_P)) + (Per_N2*(xN2*P_F + yN2*P_P)) + (Per_O2*(xO2*P_F + yO2*P_P)) );

eq2 = (xH2- xH2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_H2*(xH2*P_F + yH2*P_P) + xH2*(Ff-Ffprev)/h);
eq3 = (xN2- xN2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_N2*(xN2*P_F + yN2*P_P) + xN2*(Ff-Ffprev)/h);
eq4 = (xNH3- xNH3prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_NH3*(xNH3*P_F + yNH3*P_P) + xNH3*(Ff-Ffprev)/h);
eq5 = (xO2- xO2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_O2*(xO2*P_F + yO2*P_P) + xO2*(Ff-Ffprev)/h);

eq6 = yH2 – (Per_H2*xH2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_H2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq7 = yN2 – (Per_N2*xN2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_N2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq8 = yNH3 – (Per_NH3*xNH3 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_NH3*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq9 = yO2 – (Per_O2*xO2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_O2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));

%eq10 = ((Pp – Ppprev)/h) -( (8*R*T*mu*(F_feed-Ff))/ (3.14*r_fiber^4*n_fiber*Pp));

F = [eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9];
J=jacobian([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9],variables);
x_0 = [10 ;0.5; 0.001; 0.1; 0.0001; 0.3; 0.01; 1; 0.001];

Final_results(1,1:9,1) = x_0′;
iterations = 100;
for iter=1:iterations
%% Newton
% Substitute the initial values into the Jacobian matrix
J_subs = subs(J, variables, x_0);
F_subs = subs(F,variables,x_0);
aug_matrix = [J_subs, -F_subs];
n= 9;
for i = 1:n-1
for j = i+1:n
factor = aug_matrix(j, i) / aug_matrix(i, i);
aug_matrix(j, 🙂 = aug_matrix(j, 🙂 – factor * aug_matrix(i, :);
end
end

y_0 = zeros(n, 1);
y_0(n) = aug_matrix(n, n+1) / aug_matrix(n,n);
for i = n-1:-1:1
y_0(i) = (aug_matrix(i, n+1) – aug_matrix(i, i+1:n) * y_0(i+1:n)) / aug_matrix(i, i);
end

x_result = y_0 + x_0;
Final_results(iter+1,1:9,z) = x_result’;
err = norm(x_result – x_0);
Final_results(iter+1,10,z) = err;
x_0 = x_result;
if err <= 1e-8

%P_Pprev = x_result(10);

disp([‘Converged after ‘, num2str(iter), ‘ iterations’]);
break;
end
end
xH2prev = x_0(2);
xN2prev = x_0(3);
xNH3prev = x_0(4);
xO2prev = x_0(5);
Ffprev = x_0(1);

end

solution_Ff_Per_step = [];
solution_xH2_Per_step = [];
solution_xN2_Per_step = [];
solution_xNH3_Per_step = [];
solution_xO2_Per_step = [];
solution_yH2_Per_step = [];
solution_yN2_Per_step = [];
solution_yNH3_Per_step = [];
solution_yO2_Per_step = [];
solution_P_P_Per_step = [];

for m=1:numel(mesh)
solution_Ff_Per_step = [Final_results(end,1,m) solution_Ff_Per_step];
solution_xH2_Per_step = [Final_results(end,2,m) solution_xH2_Per_step];
solution_xN2_Per_step = [Final_results(end,3,m) solution_xN2_Per_step];
solution_xNH3_Per_step = [Final_results(end,4,m) solution_xNH3_Per_step];
solution_xO2_Per_step = [Final_results(end,5,m) solution_xO2_Per_step];
solution_yH2_Per_step = [Final_results(end,6,m) solution_yH2_Per_step];
solution_yN2_Per_step = [Final_results(end,7,m) solution_yN2_Per_step];
solution_yNH3_Per_step = [Final_results(end,8,m) solution_yNH3_Per_step];
solution_yO2_Per_step = [Final_results(end,9,m) solution_yO2_Per_step];
%solution_P_P_Per_step = [Final_results(end,10,m) solution_P_P_Per_step];

end
figure(1)
plot(mesh,solution_Ff_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Flow rate retentate [m/s]’)
figure(2)
plot(mesh,solution_xH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in retentate’)
figure(3)
plot(mesh,solution_xN2_Per_step )
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in retentate’)
figure(4)
plot(mesh,solution_xNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in retentate’)
figure(5)
plot(mesh,solution_xO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in retentate’)
figure(6)
plot(mesh,solution_yH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in Permeate’)
figure(7)
plot(mesh,solution_yN2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in Permeate’)
figure(8)
plot(mesh,solution_yNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in Permeate’)
figure(9)
plot(mesh,solution_yO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in Permeate’)
figure(10)
% plot(mesh,solution_P_P_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Permeate pressure [bar]’)
For some reason, when the method converges it will always return the initial values of xH2prev, xNH3prev, xO2prev. xN2prev and Ffprev, meaning that there is no change in mole fractions along the membrane. When I try to change the initial conditions, the graphs show spikes in random positions whilst being 0 in other positions along the membrane. This should not be the case. Could you please help me figure out what went wrong in my code? iterative, newton, backwards difference, finite MATLAB Answers — New Questions

​

Hello everyone, I would like to use pdepe for solving a SIR partial differential equations in one dimension with time lag, so it looks good. But there has been a problem setting the value of u_tau, resulting in a very different result than expected. I don’t know what’s wrong with the code.Here is the part of the code involved：
function u_tau = interpolate_history(history, x, t, tau)
% Interpolating function to get the state at time t – tau
t_target = t – tau;
if t_target < history.t(1)
% u_tau = squeeze(history.u(1, :, :)); % Out of history range, use first value
u_tau = 0;
else
% u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, ‘linear’, ‘extrap’);
u_tau = 0;
end
end

function u0 = SIRInitialConditions3(x)
S0 = 1.25;
I0 = 0.85;
R0 = 0.74;
u0 = [S0; I0; R0];
end

% This is the part of the main function that deals with u_tau
function [c,f,s] = SIRPDE3(x, t, u, DuDx)
global history;
S = u(1);
I = u(2);
R = u(3);
u_tau = interpolate_history(history, x, t, tau);
S_pre = u_tau(1);
end

% This is the part of the main function that deals with u_tau
function run()
xspan = linspace(0, 10, 100);
tspan = [0, 10, 50];

% Initialising the history
history = initialize_history(xspan, tspan, @SIRInitialConditions3);

sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
I originally ran the simulation with the two equations commented above neither of which resulted in a trend over time, so I set all u_tau to 0. The trend over time then turned out to be correct, but the initial values of I and R were still 0, which was not the expected result. I think there might be something wrong with the part function u_tau = interpolate_history(history, x, t, tau), how should I change it?
Thanks to the community. :)Hello everyone, I would like to use pdepe for solving a SIR partial differential equations in one dimension with time lag, so it looks good. But there has been a problem setting the value of u_tau, resulting in a very different result than expected. I don’t know what’s wrong with the code.Here is the part of the code involved：
function u_tau = interpolate_history(history, x, t, tau)
% Interpolating function to get the state at time t – tau
t_target = t – tau;
if t_target < history.t(1)
% u_tau = squeeze(history.u(1, :, :)); % Out of history range, use first value
u_tau = 0;
else
% u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, ‘linear’, ‘extrap’);
u_tau = 0;
end
end

function u0 = SIRInitialConditions3(x)
S0 = 1.25;
I0 = 0.85;
R0 = 0.74;
u0 = [S0; I0; R0];
end

% This is the part of the main function that deals with u_tau
function [c,f,s] = SIRPDE3(x, t, u, DuDx)
global history;
S = u(1);
I = u(2);
R = u(3);
u_tau = interpolate_history(history, x, t, tau);
S_pre = u_tau(1);
end

% This is the part of the main function that deals with u_tau
function run()
xspan = linspace(0, 10, 100);
tspan = [0, 10, 50];

% Initialising the history
history = initialize_history(xspan, tspan, @SIRInitialConditions3);

sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
I originally ran the simulation with the two equations commented above neither of which resulted in a trend over time, so I set all u_tau to 0. The trend over time then turned out to be correct, but the initial values of I and R were still 0, which was not the expected result. I think there might be something wrong with the part function u_tau = interpolate_history(history, x, t, tau), how should I change it?
Thanks to the community. 🙂 Hello everyone, I would like to use pdepe for solving a SIR partial differential equations in one dimension with time lag, so it looks good. But there has been a problem setting the value of u_tau, resulting in a very different result than expected. I don’t know what’s wrong with the code.Here is the part of the code involved：
function u_tau = interpolate_history(history, x, t, tau)
% Interpolating function to get the state at time t – tau
t_target = t – tau;
if t_target < history.t(1)
% u_tau = squeeze(history.u(1, :, :)); % Out of history range, use first value
u_tau = 0;
else
% u_tau = interp1(history.t, squeeze(history.u(:, :, :)), t_target, ‘linear’, ‘extrap’);
u_tau = 0;
end
end

function u0 = SIRInitialConditions3(x)
S0 = 1.25;
I0 = 0.85;
R0 = 0.74;
u0 = [S0; I0; R0];
end

% This is the part of the main function that deals with u_tau
function [c,f,s] = SIRPDE3(x, t, u, DuDx)
global history;
S = u(1);
I = u(2);
R = u(3);
u_tau = interpolate_history(history, x, t, tau);
S_pre = u_tau(1);
end

% This is the part of the main function that deals with u_tau
function run()
xspan = linspace(0, 10, 100);
tspan = [0, 10, 50];

% Initialising the history
history = initialize_history(xspan, tspan, @SIRInitialConditions3);

sol = pdepe(0, @SIRPDE3, @SIRInitialConditions3, @SIRBoundaryConditions3, xspan, tspan);
I originally ran the simulation with the two equations commented above neither of which resulted in a trend over time, so I set all u_tau to 0. The trend over time then turned out to be correct, but the initial values of I and R were still 0, which was not the expected result. I think there might be something wrong with the part function u_tau = interpolate_history(history, x, t, tau), how should I change it?
Thanks to the community. 🙂 pdepe, lag, diffusion equation, sir model MATLAB Answers — New Questions

​

Hi,
I am testing my environment using the PPO algorithm in RL toolbox, I recently viewed this paper: https://arxiv.org/abs/2012.01399 which listed some assumptions on the convergence guranteen of PPO, some of them are for the environment itself (like the transition kernel…) and some are for the functions and parameters of the algorithm (like the learning rate alpha, the update function h…)
I am not sure if the PPO algorithm in the RL toolbox satisfies the assumptions of the convergence for the functions and parameters of the algorithm, because I did not find any direct mentioning of convergence in the official mathwork website, so I wonder how the algorithm is designed such that convergence is being considered.
Do I need to look into the train() function to see how those parameters and functions are designed?
Thank youHi,
I am testing my environment using the PPO algorithm in RL toolbox, I recently viewed this paper: https://arxiv.org/abs/2012.01399 which listed some assumptions on the convergence guranteen of PPO, some of them are for the environment itself (like the transition kernel…) and some are for the functions and parameters of the algorithm (like the learning rate alpha, the update function h…)
I am not sure if the PPO algorithm in the RL toolbox satisfies the assumptions of the convergence for the functions and parameters of the algorithm, because I did not find any direct mentioning of convergence in the official mathwork website, so I wonder how the algorithm is designed such that convergence is being considered.
Do I need to look into the train() function to see how those parameters and functions are designed?
Thank you Hi,
I am testing my environment using the PPO algorithm in RL toolbox, I recently viewed this paper: https://arxiv.org/abs/2012.01399 which listed some assumptions on the convergence guranteen of PPO, some of them are for the environment itself (like the transition kernel…) and some are for the functions and parameters of the algorithm (like the learning rate alpha, the update function h…)
I am not sure if the PPO algorithm in the RL toolbox satisfies the assumptions of the convergence for the functions and parameters of the algorithm, because I did not find any direct mentioning of convergence in the official mathwork website, so I wonder how the algorithm is designed such that convergence is being considered.
Do I need to look into the train() function to see how those parameters and functions are designed?
Thank you reinforcement learning, ppo, convergence MATLAB Answers — New Questions

​

hi…can someone help me with this question?? i need to calculate the summation of second column in both arrays…thank uhi…can someone help me with this question?? i need to calculate the summation of second column in both arrays…thank u hi…can someone help me with this question?? i need to calculate the summation of second column in both arrays…thank u summation, arrays MATLAB Answers — New Questions

​

How can I connect the upper endpoint of the red curve to the red dot while keeping the curvature nearly the same?

clear all
format long
warning off
figure
set(0,’DefaultAxesFontSize’,20);
load(‘BT_Hom(2).mat’)
[~,idx] = max(x(326,:));
plot(x(326,1:idx-20),x(325,1:idx-20),’r’, ‘LineWidth’,2);
hold on
axis([0.14 .18 .15 .18]);
scatter(0.174701,0.177614, ‘o’, ‘MarkerFaceColor’, ‘r’);
xlabel(‘$arightarrow$’,’FontSize’,20,’interpreter’,’latex’,’FontWeight’,’normal’,’Color’,’k’)
ylabel(‘$brightarrow$’,’FontSize’,20,’interpreter’,’latex’,’FontWeight’,’normal’,’Color’,’k’)How can I connect the upper endpoint of the red curve to the red dot while keeping the curvature nearly the same?

clear all
format long
warning off
figure
set(0,’DefaultAxesFontSize’,20);
load(‘BT_Hom(2).mat’)
[~,idx] = max(x(326,:));
plot(x(326,1:idx-20),x(325,1:idx-20),’r’, ‘LineWidth’,2);
hold on
axis([0.14 .18 .15 .18]);
scatter(0.174701,0.177614, ‘o’, ‘MarkerFaceColor’, ‘r’);
xlabel(‘$arightarrow$’,’FontSize’,20,’interpreter’,’latex’,’FontWeight’,’normal’,’Color’,’k’)
ylabel(‘$brightarrow$’,’FontSize’,20,’interpreter’,’latex’,’FontWeight’,’normal’,’Color’,’k’) How can I connect the upper endpoint of the red curve to the red dot while keeping the curvature nearly the same?

clear all
format long
warning off
figure
set(0,’DefaultAxesFontSize’,20);
load(‘BT_Hom(2).mat’)
[~,idx] = max(x(326,:));
plot(x(326,1:idx-20),x(325,1:idx-20),’r’, ‘LineWidth’,2);
hold on
axis([0.14 .18 .15 .18]);
scatter(0.174701,0.177614, ‘o’, ‘MarkerFaceColor’, ‘r’);
xlabel(‘$arightarrow$’,’FontSize’,20,’interpreter’,’latex’,’FontWeight’,’normal’,’Color’,’k’)
ylabel(‘$brightarrow$’,’FontSize’,20,’interpreter’,’latex’,’FontWeight’,’normal’,’Color’,’k’) plot, fiigure MATLAB Answers — New Questions

​