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

조회 수: 5(최근 30일)
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
Hashim
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.

댓글을 달려면 로그인하십시오.

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 3월 17일
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
Mathieu NOE
Mathieu NOE 2021년 4월 6일
you can do the same without a for loop
ind = find(abs(ddy(i)) < 1e-9);
disp(t(ind))

댓글을 달려면 로그인하십시오.

추가 답변(0개)

범주

Find more on Delaunay Triangulation 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!

Translated by