How to find the steady state using events function for pde?

3 次查看（过去 30 天）

显示更早的评论

Hashim 2025-2-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2174389-how-to-find-the-steady-state-using-events-function-for-pde

编辑： Torsten 2025-2-23

在 MATLAB Online 中打开

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Hashim 2025-2-23

Why is that? Can you not take percent error between two subsequent readings for dc/dt or d2c/dx2?

Then we can look at the delta of the time change or spatial change. Surely it is a flawed approach as error oscillations come in for numerical solutions but I do not see why you would say that it is not useful. What makes it useful then if not this?

Unfortunatley it is quite hard to solve for unsteady PDEs otherwise we would have time constant for PDEs?

What would be an ideal approach according to you?

Torsten 2025-2-23

编辑：Torsten 2025-2-23

You are given the solution at time t and the corresponding mesh values in the event function. Can you can tell me a senseful indicator from these values that steady-state has reached ?

As ODE15S approaches steady-state, the time steps it makes become very large. Assume steady-state is reached at t = 40. Then whether you specify tend = 50 or tend = 50000 will not make a difference in computation time. So why make a fuss of it ?

请先登录，再进行评论。

请先登录，再回答此问题。