Comparing two row vectors to find constant slope for steady state condition
4 views (last 30 days)
Show older comments
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 Comments
Mathieu NOE
on 17 Mar 2021
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 ?
Accepted Answer
Mathieu NOE
on 17 Mar 2021
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/554127/image.png)
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 Comments
Mathieu NOE
on 6 Apr 2021
you can do the same without a for loop
ind = find(abs(ddy(i)) < 1e-9);
disp(t(ind))
More Answers (0)
See Also
Categories
Find more on Audio Processing Algorithm Design in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!