How to get the time step used in PDEPE solver?

For PDEPE, Matlab guide says. "The pdepe function performs the time integration with an ODE solver that selects both the time step and the formula dynamically."

How to get the actual time step used by PDEPE as an ouput?

Answer 1

[Torsten]

pdepe uses ode15s as solver for the arising differential equations.

This solver uses adaptive time stepping - thus if the problem is difficult, small time steps are chosen, if the problem is easy, large time steps are chosen.

The basic time steps (time steps that led to a successful integration step because the error tolerance was satisfied) can be seen after integration has finished if you choose only start and end point of integration for the definition of tspan:

tspan = [tstart tend]

Getting time steps to modify your equations during integration is almost impossible.

Could you explain why you want to get the time steps as output ?

[Gopinath NK]

@Torsten Thanks for clarifying. I am not planning to modify the equations during integration.

The reason for wanting to get the time steps as output is:

• To see variation in certain field quantities (like heat flux) with time. While I can get such field quantities as output, the array size for such field quantities is not consistent with the user-defined tspan. So I need to know the time steps used by ode15s solver in order to calculate the physical time marching.

Attached is a plot showing temperature vs. heat flux. As it can be seen, the ode15s is choosing variable sizes of time steps (by looking at the cluster of points). Essentially I am looking to plot time (tspan) vs. heat flux.

Is it possible?

[Torsten]

But as far as I know, heat flux is lambda*A*dT/dx.

What does it have to do with the time step ?

What you get from pdepe is a time vector and a vector of temperatures over the domain.

From the temperatures, you can calculate heat flux from the above formula for all times and plot it against temperature.

[Gopinath NK]

lambda*A*dT/dx is just the conduction part of the heat flux in the governing heat equation.

At the domain boundaries as we specify convective or radiative flux, the resultant heat flux in the system will differ.

I am interested in knowing this 'resultant' heat flux variation with time as the solution marches in time and space, which cannot be lambda*A*dT/dx.

Anyway, getting back to the original question, can we probe the variable time steps used by ode15s and save it as ouput? This output would solve my question.

PS: I found another similar question, but it was not answered (https://in.mathworks.com/matlabcentral/answers/434750-how-to-find-the-exact-integration-timesteps-of-the-pdepe-function)

[Torsten]

None of the heat fluxes you list is explicitly dependent on time - all can be calculated directly from the temperatures. If inputs to radiative or convective heat flux depend on time (external temperature, heat transfer coefficients,...), it's of course different.

Anyway:

You call pdepe with a tspan vector, and the solution pdepe returns is the solution at the time instants you have specified.

After looking closer at pdepe, access to the internal time stepping seems only possible by an OutputFcn function that you define. Since this is different to ODE15S, you should check first with a simple test code whether OutputFcn is really called after each successful time step or also only at the time instants that you specified in tspan. In the code below, function out is not called from pdepe - I don't know why.

x = linspace(0,1,25);
t = linspace(0,1,25);
m = 1;
options = odeset('OutputFcn',@out);
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t,options); 

function status = out(t,y,flag)
disp(flag)
status = 0;
end
function [c,f,s] = heatcyl(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end

[Gopinath NK]

Thanks. I shall give this a try.

"If inputs to radiative or convective heat flux depend on time (external temperature, heat transfer coefficients,...), it's of course different" - Yes, this is my case.

[Torsten]

I'd simply use "tspan" according to your needs and calculate the heat fluxes from the solution "sol" (which is supplied for the times defined in "tspan"). Then you don't need the "OutputFcn".