How to find the steady state using events function for pde?
Hi,
I am solving coupled reaction diffusion PDEs of the form:
I am using the ‘pdepe’ function to solve them. I have become aware that there is an ‘events’ functionality embedded into the ‘pdepe’ solver that allows the user to trigger an event when a certain event occurs. This is because the pdepe solver uses the ODE15s solver for dynamic time integration. I would like to use this functionality to find when the ‘steady state’ event occurs. I believe there are limited examples when using this functionality for PDEs.
I am using the following code to call ‘pdepe’ in the function file:
optns = odeset(‘Events’,@ssEvent);
[sol1,tsol,sole,te,ie] = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, alpha,…
kappa, eta, gamma, mu),IC, BC, chi, t, optns);
I then resolved the ‘sol1’ into following components.
% Concentration Profiles c(i, j, k)(Solutions)
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
Standard syntax for PDEs for event function is as:
function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
One of my questions is, is the ‘umesh’ the same as ‘sol1’?
Do I have to provide a matrix to ‘isterminal’ to inform the function to terminate or just a single value is enough? Same for direction.
Currently I am using the following script to try to find the steady state. But it is not working very well so I am here.
function [value, isterminal, direction] = ssEvent(m, t, xmesh, umesh)
% Compute the absolute difference between successive time steps
diff_solution = abs(diff(umesh, 1, 2)); % Max change in solution across time
% Define steady-state condition (threshold for change)
steady_threshold = 1e-4;
% Event occurs when the maximum solution change is below the threshold
value = diff_solution – steady_threshold;
isterminal = 0; % Stop integration when steady state is reached
direction = 0; % Detects steady state in both increasing & decreasing directions
end
Eventual goal is to make a graph like this. Where I want to extract the time against the event and plot it against the gamma. This graph was made using another method though.
Also I would appreciate if respected members can give some pointers on how to rid of numerical oscillations while computing errors.Hi,
I am solving coupled reaction diffusion PDEs of the form:
I am using the ‘pdepe’ function to solve them. I have become aware that there is an ‘events’ functionality embedded into the ‘pdepe’ solver that allows the user to trigger an event when a certain event occurs. This is because the pdepe solver uses the ODE15s solver for dynamic time integration. I would like to use this functionality to find when the ‘steady state’ event occurs. I believe there are limited examples when using this functionality for PDEs.
I am using the following code to call ‘pdepe’ in the function file:
optns = odeset(‘Events’,@ssEvent);
[sol1,tsol,sole,te,ie] = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, alpha,…
kappa, eta, gamma, mu),IC, BC, chi, t, optns);
I then resolved the ‘sol1’ into following components.
% Concentration Profiles c(i, j, k)(Solutions)
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
Standard syntax for PDEs for event function is as:
function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
One of my questions is, is the ‘umesh’ the same as ‘sol1’?
Do I have to provide a matrix to ‘isterminal’ to inform the function to terminate or just a single value is enough? Same for direction.
Currently I am using the following script to try to find the steady state. But it is not working very well so I am here.
function [value, isterminal, direction] = ssEvent(m, t, xmesh, umesh)
% Compute the absolute difference between successive time steps
diff_solution = abs(diff(umesh, 1, 2)); % Max change in solution across time
% Define steady-state condition (threshold for change)
steady_threshold = 1e-4;
% Event occurs when the maximum solution change is below the threshold
value = diff_solution – steady_threshold;
isterminal = 0; % Stop integration when steady state is reached
direction = 0; % Detects steady state in both increasing & decreasing directions
end
Eventual goal is to make a graph like this. Where I want to extract the time against the event and plot it against the gamma. This graph was made using another method though.
Also I would appreciate if respected members can give some pointers on how to rid of numerical oscillations while computing errors. Hi,
I am solving coupled reaction diffusion PDEs of the form:
I am using the ‘pdepe’ function to solve them. I have become aware that there is an ‘events’ functionality embedded into the ‘pdepe’ solver that allows the user to trigger an event when a certain event occurs. This is because the pdepe solver uses the ODE15s solver for dynamic time integration. I would like to use this functionality to find when the ‘steady state’ event occurs. I believe there are limited examples when using this functionality for PDEs.
I am using the following code to call ‘pdepe’ in the function file:
optns = odeset(‘Events’,@ssEvent);
[sol1,tsol,sole,te,ie] = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, alpha,…
kappa, eta, gamma, mu),IC, BC, chi, t, optns);
I then resolved the ‘sol1’ into following components.
% Concentration Profiles c(i, j, k)(Solutions)
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
Standard syntax for PDEs for event function is as:
function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
One of my questions is, is the ‘umesh’ the same as ‘sol1’?
Do I have to provide a matrix to ‘isterminal’ to inform the function to terminate or just a single value is enough? Same for direction.
Currently I am using the following script to try to find the steady state. But it is not working very well so I am here.
function [value, isterminal, direction] = ssEvent(m, t, xmesh, umesh)
% Compute the absolute difference between successive time steps
diff_solution = abs(diff(umesh, 1, 2)); % Max change in solution across time
% Define steady-state condition (threshold for change)
steady_threshold = 1e-4;
% Event occurs when the maximum solution change is below the threshold
value = diff_solution – steady_threshold;
isterminal = 0; % Stop integration when steady state is reached
direction = 0; % Detects steady state in both increasing & decreasing directions
end
Eventual goal is to make a graph like this. Where I want to extract the time against the event and plot it against the gamma. This graph was made using another method though.
Also I would appreciate if respected members can give some pointers on how to rid of numerical oscillations while computing errors. pdepe, ode15s MATLAB Answers — New Questions