%% Solve Oscillatory PDE with Event Logging
% Solve a partial differential equation and use an event function to log zero-crossings in the oscillatory solution.
% Consider the equation
% \frac{1}{x}\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(\frac{1}{t}u\right)\ldotp
% The equation is defined for 0\le x\le 1 and t\ge 0\ldotp 1. The initial condition is
% u\left(x,0\ldotp 1\right)=1\ldotp
% The boundary conditions are
% u\left(0,t\right)=1,
% u\left(1,t\right)=\cos \left(\pi t\right)\ldotp
% Additionally, the zero-crossings of the solution are of interest.
% To solve this equation in MATLAB, you need to code the equation, initial conditions, boundary conditions, and event function, then select a suitable solution mesh before calling the solver pdepe. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.
%% Code Equation
% Before you can code the equation, you need to rewrite it in a form that the pdepe solver expects. The standard form that pdepe expects is
% c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}=x^{-m} \frac{\partial }{\partial x}\left(x^m f\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right)\ldotp
% The PDE equation is already in this form:
% \frac{1}{x}\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(\frac{1}{t}u\right)\ldotp
% So you can read off the relevant terms:
% m=0
% c\left(x,t,u,\frac{\partial u}{\partial x}\right)=\frac{1}{x}
% f\left(x,t,u,\frac{\partial u}{\partial x}\right)=\frac{1}{t}u
% s\left(x,t,u,\frac{\partial u}{\partial x}\right)=0
% Now you can create a function to code the equation. The function should have the signature [c,f,s] = oscpde(x,t,u,dudx):
% x is the independent spatial variable.
% t is the independent time variable.
% u is the dependent variable being differentiated with respect to x and t.
% dudx is the partial spatial derivative \partial u/\partial x.
% The outputs c, f, and s correspond to coefficients in the standard PDE equation form expected by pdepe. These coefficients are coded in terms of the input variables x, t, u, and dudx.
% As a result, the equation in this example can be represented by the function:
% function [c,f,s] = oscpde(x,t,u,dudx)
% c = 1/x;
% f = u/t;
% s = 0;
% end
% (Note: All functions are included as local functions at the end of the example.)
%
%% Code Initial Conditions
% Next, write a function that returns the initial condition. The initial condition is applied at the first time value tspan(1). The function should have the signature u0 = oscic(x).
% The corresponding function for u\left(x,0\ldotp 1\right)=1 is
% function u0 = oscic(x)
% u0 = 1;
% end
%
%% Code Boundary Conditions
% Now, write a function that evaluates the boundary conditions
% u\left(0,t\right)=1,
% u\left(1,t\right)=\cos \left(\pi t\right)\ldotp
% For problems posed on the interval a\le x\le b, the boundary conditions apply for all t and either x=a or x=b. The standard form for the boundary conditions expected by the solver is
% p\left(x,t,u,\frac{\partial u}{\partial x}\right)+q\left(x,t\right)f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0\ldotp
% Written in this form, the boundary conditions for this problem become
% u\left(0,t\right)-1+0\cdot f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0,
% u\left(1,t\right)-\cos \left(\pi x\right)+0\cdot f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0\ldotp
% The boundary function should have the function signature [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t):
% The inputs xl and ul correspond to u and x for the left boundary.
% The inputs xr and ur correspond to u and x for the right boundary.
% t is the independent time variable.
% The outputs pl and ql correspond to p_L \left(x,t,u\right) and q_L \left(x,t\right) for the left boundary (x=0 for this problem).
% The outputs pr and qr correspond to p_R \left(x,t,u\right) and q_R \left(x,t\right) for the right boundary (x=1 for this problem).
%  The boundary conditions in this example are represented by the function:
% function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
% pl = ul - 1;
% ql = 0;
% pr = ur - cos(pi*t);
% qr = 0;
% end
%
%% Code Event Function
% Use an event function to log the zero-crossings of the solution in the integration. The event function has a function signature of [value,isterminal,direction] = pdevents(m,t,xmesh,umesh):
% m is the coordinate symmetry specified as the first input to pdepe.
% t is the current time (a scalar).
% xmesh is the spatial mesh.
% umesh contains the solution at the mesh points.
% value is an equation of interest, usually expressed in terms of the solution umesh. When value equals 0, an event occurs.
% isterminal specifies whether an event stops the integration. If isterminal is 0, events are logged but the integration does not stop. If isterminal is 1, then the integration stops when an event occurs.
% direction specifies the direction of zero crossing. If 1, only zero-crossings with a positive slope trigger an event. If -1, the zero-crossings must have a negative slope. If 0, then any zero-crossing triggers an event.
% At each time in the integration, the solver calls the event function to check for zero crossings. To log all zero crossings, value should look for changes in sign in the solution vector umesh. Specify isterminal and direction as vectors of zeros with the same size, since the events in this example are not terminal and the zero-crossings can occur with any slope.
% The event function for this problem is
% function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
% value = umesh;
% isterminal = zeros(size(umesh));
% direction = zeros(size(umesh));
% end
%
%% Select Solution Mesh
% Before solving the equation you need to specify the mesh points \left(t,x\right) at which you want pdepe to evaluate the solution. For this problem, use a fine mesh of 50 points in the intervals 0\le x\le 1 and 0\ldotp 1\le t\le 1. A fine mesh gives good resolution of the oscillatory solution.

x = linspace(0,1,50);
t = linspace(0.1,pi,50);

%
%% Solve Equation
% Finally, solve the equation using the symmetry m, the PDE equation, the initial conditions, the boundary conditions, the event function, and the meshes for x and t. Use odeset to create an options structure that references the events function, and pass in the  structure as the last input argument to pdepe. Specify five output arguments to return information from both the event function and the solver:
% sol is the solution computed by pdepe.
% tsol is a vector of times before a terminal event. tsol is equal to t when no events are terminal.
% sole is the solution at the time of each event.
% te is the time of each event.
% ie is the index of each event. Since values = umesh in the event function, ie gives the indices of umesh that triggered an event at each time step.

m = 0;
options = odeset('Events',@pdevents);
[sol,tsol,sole,te,ie] = pdepe(m,@oscpde,@oscic,@oscbc,x,t,options);

% Extract the solution as a matrix u.

u = sol(:,:,1);

%
%% Plot Solution
% Create a surface plot of the solution and view the plot from above.

surf(x,t,u)
view(2)

%
% Plot the points where events occurred, with a surface f\left(x,t\right)=0 for reference. The output index vector ie is useful to pick out the event locations. The expression x(ie)' gives the x-values where events occurred, and the expression sole(x==x(ie)') gives the corresponding solution values.

view([39 30])
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
hold on
plot3(x(ie)',te,sole(x==x(ie)'),'r*')
surf(x,t,zeros(size(u)),'EdgeColor','flat')
hold off

%
%% Local Functions
% Listed here are the local helper functions that the PDE solver pdepe calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function [c,f,s] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end
%----------------------------------------------
function u0 = oscic(x)
u0 = 1;
end
%----------------------------------------------
function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end
%----------------------------------------------
function [value, isterminal, direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
%----------------------------------------------

% Copyright 2019 The MathWorks, Inc.