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

조회 수: 5(최근 30일)
Hashim 2021년 3월 17일
댓글: Mathieu NOE 2021년 4월 6일
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표시숨기기 이전 댓글 수: 1
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 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표시숨기기 이전 댓글 수: 11
Mathieu NOE 2021년 4월 6일
you can do the same without a for loop
ind = find(abs(ddy(i)) < 1e-9);
disp(t(ind))

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

R2017b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!