Comparing two row vectors to find constant slope for steady state condition
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hashim
am 17 Mär. 2021
Kommentiert: Mathieu NOE
am 6 Apr. 2021
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 Kommentare
Mathieu NOE
am 17 Mär. 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 ?
Akzeptierte Antwort
Mathieu NOE
am 17 Mär. 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 Kommentare
Mathieu NOE
am 6 Apr. 2021
you can do the same without a for loop
ind = find(abs(ddy(i)) < 1e-9);
disp(t(ind))
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Audio Processing Algorithm Design finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!