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.

[Torsten]

You get the solution "umesh" at time t in the grid points "xmesh" as input argument to the events function. So you can't compute a "difference between successive time steps".

And even if the results of two subsequent time steps were available: you cannot set an absolute threshold. If the last time step was very short, results won't change much anyway. So a small difference wouldn't indicate that steady-state is reached.

[Hashim]

There are are two components to my solution as far as 'sol1' is concerned. Do I need to take that into account? Is it possible to view 'umesh' somehow?

I suppose the idea here is that is that if:

value = umesh;

It will trigger an event wherever the umesh = 0. Also umesh is itself the time component so it is not a resolved matrix like c1 and c2 from 'sol1'?

[Torsten]

umesh will be of size 2 x (number of mesh points) or (number of mesh points) x 2 - I don't know. You can simply add the line

size(umesh)

in the events function to get this information.

If you use

value = umesh;
isterminal = ones(size(umesh));
direction = zeros(size(umesh));

the solver will stop if at least one of the solution components passes 0 in at least one mesh point. But I don't understand the relation to a steady-state condition.

[Hashim]

Is it because the events function will merely raise a flag at the point when the solution umesh passes through zero. But if I set:

value = umesh(i+1)-umesh(i);

Then it does not work because umesh is the solution on the spatial grid and not the termporal array.

Also I cannot do this along the temporal dimesnion:

value = c1(i+1)-c1(i);

Because I only get the c1 after pdepe has executed?

[Torsten]

As said: I don't see a useful criterion that could indicate whether you have reached steady-state even if it were possible to access the solutions at two subsequent times.

[Hashim]

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]

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 ?