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

조회 수: 2 (최근 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)
  댓글 수: 2
Mathieu NOE
Mathieu NOE 2021년 3월 17일
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.

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

채택된 답변

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]
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))/...
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);
xl = [l-x_k, l-x_k];
yl = [-M_bulk*2, M_bulk*2];
line(xl, yl, 'Color','g','LineStyle','--','LineWidth',2);
hold off;
%% plot of first and second derivative of C2
[dc2, ddc2] = firstsecondderivatives(x,c2( SS,: ));
subplot(3,1,1),plot(x,c2( SS,: ))
xlabel('x (cm)'); ylabel('C2');
xlabel('x (cm)'); ylabel('dC2/dx');
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);
for i = 1:N
p(:,i) = c1(:,i)-c1(:,i+1);
if p(:,i) < 1e-10
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)] ;
%% Initial Condition
function c0 = pdev7ic(x)
global S_bulk M_bulk
c0 = [M_bulk; S_bulk; M_bulk];
%% 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;
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];
%% Analytical Portion Solution
  댓글 수: 12
Hashim 2021년 4월 6일
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
Mathieu NOE
Mathieu NOE 2021년 4월 6일
you can do the same without a for loop
ind = find(abs(ddy(i)) < 1e-9);

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

추가 답변 (0개)


Help CenterFile Exchange에서 Audio Processing Algorithm Design에 대해 자세히 알아보기




Community Treasure Hunt

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

Start Hunting!

Translated by