Comparing two row vectors to find constant slope for steady state condition

Question

Hashim 2021-3-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/775252-comparing-two-row-vectors-to-find-constant-slope-for-steady-state-condition

评论： Mathieu NOE 2021-4-6

pdepe_philip_v7a.m

I am comparing two rows to find if at any point I am getting a uniform slope (a condition for steady state let's say in a diffusion problem)

c2(N+1,N+1)=9.90083146790674e-07;    % Doing this to add an extra index so I can loop through this 100*100 matrix 
for i = 1:N 
    p(:,i) = c2(:,i)-c2(:,i+1); 
    if p(:,i) < 1e-10         % I wany to see if the slope is uniform/0 so I can ascertain if my solution has reached steady state (for a diffusion problem) 
        disp(i)
    end
end

Thanks!

2 个评论
显示无隐藏无

Mathieu NOE 2021-3-17

hello

so basically you are doing a difference (like diff) and compare that to a threshold.

That's a valid approach to find a constant slope condition as you say (steady state)

so what is the issue ? what are you willing to do ?

Hashim 2021-3-17

编辑：Hashim 2021-3-17

Hi Sorry I got the wrong variable here... it's actually c2 (see attached file) that I am concerned with. if you run the code you will find it output a graph (figure 2) with two y axis. I wish to know at what time the slope on the right curve is removed or lessened signifncatly so I can know that my system is at a steady state.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2021-3-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/775252-comparing-two-row-vectors-to-find-constant-slope-for-steady-state-condition#answer_650562

在 MATLAB Online 中打开

firstsecondderivatives.m

Hi again

so I introduced the first and second derivative of c2 and plotted in figure 3

Attached is the function to do a finite difference derivation better than using diff (or alike)

as the second derivative goes to zero in the second half of the x data, you can tell the system has reached a steady state

main code (modified) below :

% function pdepe_philip_v7a
% 5- We now add a second order reaction to the mix and try to study the
% effects of this reaction on the concentration output of the system 
% Also, to understand the concept of reaction layers and how they can explain 
% behavior or our system 
% Solve for Substrate concentration as well 
global D_M D_S M_bulk n F m Area R S_bulk k1
n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D_M        = 5e-06;                              % cm^2/s 
D_S        = 5e-06;                              % cm^2/s (???)
l          = 0.01;                               % cm 
Area       = 1;                                  % cm^2 
M_bulk     = 1e-06;                              % mol/cm^3
S_bulk     = 1e-04;                              % mol/cm^3 
m          = 0;                                  % Cartesian Co-ordinates
k1         = 1e+05; 
N          = 100; 
 
% computational cost of the solution depends weakly on the length of tspan
t          = linspace(0, 10, N);                % s       
% xmesh (Determines the computational cost)
x          = linspace(0, l, N);
sol        = pdepe(m, @pdev7, @pdev7ic, @pdev7bc, x, t);
c1         = sol(:, :, 1);                        % Mox Conc.  
c2         = sol(:, :, 2);                        % Substrate Conc.
c3         = sol(:, :, 3);                        % Mred Conc.  
% % Conc. Profiles (3D)
% figure(1)
% surf(x,t,c2);       
% view(0, 0);
% xlabel('x');
% ylabel('t(sec)');
% zlabel('Concentration');
% % dc/dt [Glucose] 
% figure(3); 
% plot(x, c1( 2:N, :));  
SS         = 94;
% dc/dt [Os] + [Glucose]
figure(2); 
plotyy(x, c1( SS, :), x, c2( SS,: )); 
xlabel('x (cm)'); ylabel('[Os^3] (mol/cm^3)'); 
hold on;  
% draw a line for sigma_kinetic/reaction layer [cm]
x_k        = sqrt(((D_S*M_bulk)+(D_M*S_bulk))/...
             (k1*(M_bulk+S_bulk).^2));     
         
if S_bulk > M_bulk
    xl          = [x_k, x_k];
    yl          = [-M_bulk*2, M_bulk*2];
    line(xl, yl, 'Color','green','LineStyle','--','LineWidth',2);
else 
    xl          = [l-x_k, l-x_k];
    yl          = [-M_bulk*2, M_bulk*2];
    line(xl, yl, 'Color','g','LineStyle','--','LineWidth',2);
end
hold off;  
%%  plot of first and second derivative of C2
  [dc2, ddc2] = firstsecondderivatives(x,c2( SS,: ));
figure(3); 
subplot(3,1,1),plot(x,c2( SS,: ))
xlabel('x (cm)'); ylabel('C2'); 
subplot(3,1,2),plot(x,dc2)
xlabel('x (cm)'); ylabel('dC2/dx'); 
subplot(3,1,3),plot(x,ddc2)
xlabel('x (cm)'); ylabel('d²C2/dx²'); 
% % Mred Flux at x=0
% figure(4);
% subplot(4,1,1);
% plot(t, c1( :,1)); 
% xlabel('time (s)'); ylabel('dM/dx @ x=0)');  
% 
% % Mred Flux at x=d
% subplot(4,1,2);
% plot(t, c1(:,N));
% xlabel('time (s)'); ylabel('dM/dx @ x=d)');  
% 
% % Substrate Flux at x=d
% subplot(4,1,3);
% plot(t, c2(:,N));
% xlabel('time (s)'); ylabel('dS/dx @ x=d)');  
% 
% % Mred Flux at x=0 - Mred Flux at x=0
% subplot(4,1,4);
% plot(t, c1(:,N)-c1(:,1));
% xlabel('time (s)'); ylabel('dM/dx@x=d-dM/dx@x=0');  
asd  = c1(:,N-1) - c1(:,N); 
c1(N+1,N+1)=9.90083146790674e-07;
for i = 1:N 
    p(:,i) = c1(:,i)-c1(:,i+1); 
    if p(:,i) < 1e-10
        disp(i)
    end
end
as = 1
% end
%% pdepe Function That Calls Children
%%
function [a, f, s] = pdev7(x, t, c, DuDx)
global D_M D_S k1  
a          = [1; 1; 1]; 
f          = [D_M; D_S; D_M].*DuDx;    
% c(1)--> Mox(x,t) || c(2)--> S(x,t) || c(3)--> Mred(x,t)
s          = [-(k1*c(1)*c(2)); -k1*c(1)*c(2); k1*c(1)*c(2)] ;
end
%% Initial Condition
%%
function c0 = pdev7ic(x)
global S_bulk M_bulk 
c0          = [M_bulk; S_bulk; M_bulk];
end
%% Boundry Conditions 
%%
function [pl, ql, pr, qr] = pdev7bc(xl, cl, xr, cr, t)
global M_bulk S_bulk n F R
E0   =0.2; 
E    =1.2; 
alpha=exp((E-E0)*((n*F)/(R*298.15)));   
pl   =[cl(1)-((M_bulk*alpha)./(1+alpha)); 0; cl(3)-(M_bulk*(1/(1+alpha)))];
ql   =[0;                                    1;                            0];
pr   =[0;                                    cr(2)-S_bulk;                 0];
qr   =[1;                                    0;                            1];   
end
%% Analytical Portion Solution 
%%

12 个评论
显示 10更早的评论隐藏 10更早的评论

Hashim 2021-4-4

编辑：Hashim 2021-4-4

在 MATLAB Online 中打开

Hi @Mathieu NOEMathew I've got the code working but I need to see the variable against time and not x axis. I've done that in the code below. Now I would like to check if it works by taking a snapshot of system when it is at an earlier time. For this I need to plot t(some earlier value of time) against c2(some earlier value of time). Now we know that c2 is a square matrix and time is a vector. How do I code it?

I am talking about these line BTW in the code below.

% plot of first and second derivative of C2

[dy, ddy] = firstsecondderivatives(t, c2(SS,:));

figure(3);

subplot(3,1,1),plot(t,c2(SS,: ));

xlabel('t (s)'); ylabel('C2');

subplot(3,1,2),plot(t,dy)

xlabel('t (s)'); ylabel('dC2/dx');

subplot(3,1,3),plot(t,ddy)

xlabel('t (s)'); ylabel('d²C2/dx²');

function pdepe_philip_v7a
% 5- We now add a second order reaction to the mix and try to study the
% effects of this reaction on the concentration output of the system 
% Also, to understand the concept of reaction layers and how they can explain 
% behavior or our system 
% Solve for Substrate concentration as well 
global D_M D_S M_bulk n F m Area R S_bulk k1
n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D_M        = 5e-06;                              % cm^2/s 
D_S        = 5e-06;                              % cm^2/s (???)
l          = 0.01;                               % cm 
Area       = 1;                                  % cm^2 
M_bulk     = 1e-06;                              % mol/cm^3
S_bulk     = 1e-04;                              % mol/cm^3 
m          = 0;                                  % Cartesian Co-ordinates
k1         = 1e+05; 
N          = 100; 
 
% computational cost of the solution depends weakly on the length of tspan
t          = linspace(0, 100, N);                % s       
% xmesh (Determines the computational cost)
x          = linspace(0, l, N);
sol        = pdepe(m, @pdev7, @pdev7ic, @pdev7bc, x, t);
c1         = sol(:, :, 1);                       % Mox Conc.  
c2         = sol(:, :, 2);                       % Substrate Conc.
c3         = sol(:, :, 3);                       % Mred Conc.  
% % Conc. Profiles (3D)
% figure(1)
% surf(x,t,c2);       
% view(0, 0);
% xlabel('x');
% ylabel('t(sec)');
% zlabel('Concentration');
SS         = 10;
% dc/dt [Os] + [Glucose]
figure(2); 
[hAx,hLine1,hLine2]= plotyy(x, c1( SS, :), x, c2( SS,: )); 
xlabel('x / cm'); 
ylabel(hAx(1),'Os^3^+ / mol cm^-^3')    % left y-axis 
ylabel(hAx(2),'Substrate / mol cm^-^3') % right y-axis
% legend('k=1e5 cm^3 mol s^-^1','D=5e-6 cm^2 s^-^1','t=100 s');
hold on;  
% draw a line for sigma_kinetic/reaction layer [cm]
x_k            = sqrt(((D_S*M_bulk)+(D_M*S_bulk))/...
                 (k1*(M_bulk+S_bulk).^2));     
         
if S_bulk > M_bulk
    x_new      = [x_k, x_k];
    y          = [-M_bulk*3, M_bulk*2];
    line(x_new, y, 'Color','green','LineStyle','--','LineWidth',2);
else 
    x_newer    = [l-x_k, l-x_k];
    y          = [-M_bulk*3, M_bulk*2];
    line(x_newer, y, 'Color','g','LineStyle','--','LineWidth',2);
end
size(c2(SS,1)) 
size(t)
% plot of first and second derivative of C2
[dy, ddy] = firstsecondderivatives(t, c2(SS,:));
figure(3); 
subplot(3,1,1),plot(t,c2(SS,: ));
xlabel('t (s)'); ylabel('C2'); 
subplot(3,1,2),plot(t,dy)
xlabel('t (s)'); ylabel('dC2/dx'); 
subplot(3,1,3),plot(t,ddy)
xlabel('t (s)'); ylabel('d²C2/dx²'); 
% % dMdx at x=0
% figure(4);
% subplot(4,1,1);
% plot(t, c1( :,1)); 
% xlabel('time (s)'); ylabel('dM/dx @ x=0)');  
% 
% % dMdx at x=d
% subplot(4,1,2);
% plot(t, c1(:,N));
% xlabel('time (s)'); ylabel('dM/dx @ x=d)');  
% 
% % dSdx at x=d
% subplot(4,1,3);
% plot(t, c2(:,N));
% xlabel('time (s)'); ylabel('dS/dx @ x=d)');  
% 
% % dMdx at x=0 - dMdx at x=0
% subplot(4,1,4);
% plot(t, c2(N,:)-c1(2,:));
% xlabel('time (s)'); ylabel('dM/dx@x=d-dM/dx@x=0');  
% asd  = c2(N-1,:) - c2(N,:); 
% 
% c2(N+1,N+1)=9.90083146790674e-07;
% 
% for i = 1:N 
%     p(i,:) = c2(i,:)-c2(i+1,:); 
%     if p(i,:) < 1e-8
%         disp(i)
%     end
% end
% slope   = gradient(c2(SS,:), linspace(0, l, N)); 
end
%% pdepe Function
%%
function [a, f, s] = pdev7(x, t, c, DuDx)
global D_M D_S k1  
a          = [1; 1; 1]; 
f          = [D_M; D_S; D_M].*DuDx;    
% c(1)--> Mox(x,t) || c(2)--> S(x,t) || c(3)--> Mred(x,t)
s          = [-(k1*c(1)*c(2)); -k1*c(1)*c(2); k1*c(1)*c(2)] ;
end
%% Initial Condition
%%
function c0 = pdev7ic(x)
global S_bulk M_bulk 
c0          = [M_bulk; S_bulk; M_bulk];
end
%% Boundry Conditions 
%%
function [pl, ql, pr, qr] = pdev7bc(xl, cl, xr, cr, t)
global M_bulk S_bulk n F R
E0   =0.2; 
E    =1.2; 
alpha=exp((E-E0)*((n*F)/(R*298.15)));   
pl   =[cl(1)-((M_bulk*alpha)./(1+alpha)); 0; cl(3)-(M_bulk*(1/(1+alpha)))];
ql   =[0;                                    1;                            0];
pr   =[0;                                    cr(2)-S_bulk;                 0];
qr   =[1;                                    0;                            1];   
end
%% Courtesy of Mathieu From MATLAB Website 
%%
function [dy, ddy] = firstsecondderivatives(x,y)
% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.
n      = length(x);
dy     = zeros;
ddy    = zeros;
% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1)  = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1)));       % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative
for i  = 2:n-1   
	dy(i)  = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
	ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
	dy(n)  = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
	ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end
%%
%%

Hashim 2021-4-6

在 MATLAB Online 中打开

Hi @Mathieu NOE I've been playing around with your function and pdeval which also provided dc/dx and found that even though the gradient is the same the values are off by a magnitude of 10s. Could you shine some light on this?

I've added another subplot where you can see this...

function pdepe_philip_v7a
% 5- We now add a second order reaction to the mix and try to study the
% effects of this reaction on the concentration output of the system 
% Also, to understand the concept of reaction layers and how they can explain 
% behavior or our system 
% Solve for Substrate concentration as well 
global D_M D_S M_bulk n F m Area R S_bulk k1
n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D_M        = 5e-06;                              % cm^2/s 
D_S        = 5e-06;                              % cm^2/s (???)
l          = 0.01;                               % cm 
Area       = 1;                                  % cm^2 
M_bulk     = 1e-06;                              % mol/cm^3
S_bulk     = 1e-04;                              % mol/cm^3 
m          = 0;                                  % Cartesian Co-ordinates
k1         = 1e+05; 
N          = 100; 
 
% computational cost of the solution depends weakly on the length of tspan
t          = linspace(0, 100, N);                % s       
% xmesh (Determines the computational cost)
x          = linspace(0, l, N);
sol        = pdepe(m, @pdev7, @pdev7ic, @pdev7bc, x, t);
c1         = sol(:, :, 1);                       % Mox Conc.  
c2         = sol(:, :, 2);                       % Substrate Conc.
c3         = sol(:, :, 3);                       % Mred Conc.  
SS         = 5; 
% dc/dt [Os] + [Glucose]
figure(2); 
[hAx,hLine1,hLine2]= plotyy(x, c1( SS, :), x, c2( SS,: )); 
xlabel('x / cm'); 
ylabel(hAx(1),'Os^3^+ / mol cm^-^3')    % left y-axis 
ylabel(hAx(2),'Substrate / mol cm^-^3') % right y-axis
% legend('k=1e5 cm^3 mol s^-^1','D=5e-6 cm^2 s^-^1','t=100 s');
 
% 
% % draw a line for sigma_kinetic/reaction layer [cm]
% x_k            = sqrt(((D_S*M_bulk)+(D_M*S_bulk))/...
%                  (k1*(M_bulk+S_bulk).^2));     
%          
% if S_bulk > M_bulk
%     x_new      = [x_k, x_k];
%     y          = [-M_bulk*3, M_bulk*2];
%     line(x_new, y, 'Color','green','LineStyle','--','LineWidth',2);
% else 
%     x_newer    = [l-x_k, l-x_k];
%     y          = [-M_bulk*3, M_bulk*2];
%     line(x_newer, y, 'Color','g','LineStyle','--','LineWidth',2);
% end
size(c2(1,1:SS)) 
size(t(1:SS))
% plot of first and second derivative of C2
% t1 = t(1:SS); 
c4 = c2(SS,:);
[dy, ddy] = firstsecondderivatives(t, c4);
figure(3); 
subplot(4,1,1),plot(t,c4);
xlabel('t (s)'); ylabel('C2'); 
subplot(4,1,2),plot(t,dy)
xlabel('t (s)'); ylabel('dC2/dx'); 
subplot(4,1,3),plot(t,ddy)
xlabel('t (s)'); ylabel('d²C2/dx²'); 
[cout,dcdx]= pdeval(m,x,c2(SS,:),x);
subplot(4,1,4),plot(t,dcdx)
xlabel('t (s)'); ylabel('dC2/dx');
% % dMdx at x=0
% figure(4);
% subplot(4,1,1);
% plot(t, c1( :,1)); 
% xlabel('time (s)'); ylabel('dM/dx @ x=0)');  
% 
% % dMdx at x=d
% subplot(4,1,2);
% plot(t, c1(:,N));
% xlabel('time (s)'); ylabel('dM/dx @ x=d)');  
% 
% % dSdx at x=d
% subplot(4,1,3);
% plot(t, c2(:,N));
% xlabel('time (s)'); ylabel('dS/dx @ x=d)');  
% 
% % dMdx at x=0 - dMdx at x=0
% subplot(4,1,4);
% plot(t, c2(N,:)-c1(2,:));
% xlabel('time (s)'); ylabel('dM/dx@x=d-dM/dx@x=0');  
% asd  = c2(N-1,:) - c2(N,:); 
% 
% c2(N+1,N+1)=9.90083146790674e-07;
% 
% for i = 1:N 
%     p(i,:) = c2(i,:)-c2(i+1,:); 
%     if p(i,:) < 1e-8
%         disp(i)
%     end
% end
% slope   = gradient(c2(SS,:), linspace(0, l, N)); 
end
%% pdepe Function
%%
function [a, f, s] = pdev7(x, t, c, DuDx)
global D_M D_S k1  
a          = [1; 1; 1]; 
f          = [D_M; D_S; D_M].*DuDx;    
% c(1)--> Mox(x,t) || c(2)--> S(x,t) || c(3)--> Mred(x,t)
s          = [-(k1*c(1)*c(2)); -k1*c(1)*c(2); k1*c(1)*c(2)] ;
end
%% Initial Condition
%%
function c0 = pdev7ic(x)
global S_bulk M_bulk 
c0          = [M_bulk; S_bulk; M_bulk];
end
%% Boundry Conditions 
%%
function [pl, ql, pr, qr] = pdev7bc(xl, cl, xr, cr, t)
global M_bulk S_bulk n F R
E0   =0.2; 
E    =1.2; 
alpha=exp((E-E0)*((n*F)/(R*298.15)));   
pl   =[cl(1)-((M_bulk*alpha)./(1+alpha)); 0; cl(3)-(M_bulk*(1/(1+alpha)))];
ql   =[0;                                    1;                            0];
pr   =[0;                                    cr(2)-S_bulk;                 0];
qr   =[1;                                    0;                            1];   
end
%% Courtesy of Mathieu From MATLAB Website 
%%
function [dy, ddy] = firstsecondderivatives(x,y)
% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.
n      = length(x);
dy     = zeros;
ddy    = zeros;
% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1)  = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1)));       % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative
for i  = 2:n-1   
	dy(i)  = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
	ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
	dy(n)  = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
	ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end
%%
%%

Hashim 2021-4-6

在 MATLAB Online 中打开

How about this, anyways I've accepeted the answer, thank you extremely much for your help so far.

for i = 1:length(t) 
    if abs(ddy(i)) < 1e-9
        disp(i)
    end
end

Mathieu NOE 2021-4-6

在 MATLAB Online 中打开

you can do the same without a for loop

ind = find(abs(ddy(i)) < 1e-9);
disp(t(ind))

请先登录，再进行评论。

Comparing two row vectors to find constant slope for steady state condition

2 个评论
显示无隐藏无

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Comparing two row vectors to find constant slope for steady state condition

2 个评论 显示 无隐藏 无

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

12 个评论
显示 10更早的评论隐藏 10更早的评论