%% Solve Heat Equation in Cylindrical Coordinates
% Solve the heat equation in cylindrical coordinates using pdepe, and plot the solution.
% In cylindrical coordinates with angular symmetry the heat equation is
% \frac{\partial u}{\partial t}=\frac{1}{x}\frac{\partial }{\partial x}\left(x\frac{\partial u}{\partial x}\right)\ldotp
% The equation is defined for 0\le x\le 1 at times t\ge 0. The initial condition is defined in terms of the bessel function J_0 \left(x\right) and its first zero n=2\ldotp 404825557695773 as
% u\left(x,0\right)=J_0 \left(\textrm{nx}\right)\ldotp
% Since this problem is in cylindrical coordinates (m = 1), pdepe automatically enforces the symmetry condition at x=0. The right boundary condition is
% u\left(1,t\right)=J_0 \left(n\right)e^{-n^2 t} \ldotp
% The initial and boundary conditions are selected to be consistent with the analytic solution to the problem,
% u\left(x,t\right)=J_0 \left(\textrm{nx}\right)e^{-n^2 t} \ldotp
% To solve this equation in MATLAB, you need to code the equation, initial conditions, and boundary conditions, 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
% Written in this form, the PDE becomes
% \frac{\partial u}{\partial t}=x^{-1} \frac{\partial }{\partial x}\left(x^1 \frac{\partial u}{\partial x}\right)
% With the equation in the proper form you can read off the relevant terms:
% m=1
% c\left(x,t,u,\frac{\partial u}{\partial x}\right)=1
% f\left(x,t,u,\frac{\partial u}{\partial x}\right)=\frac{\partial u}{\partial x}
% 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] = heatcyl(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] = heatcyl(x,t,u,dudx)
% c = 1;
% f = dudx;
% 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 = heatic(x).
% The corresponding function for u\left(x,0\right)=J_0 \left(\textrm{nx}\right) is
% function u0 = heatic(x)
% n = 2.404825557695773;
% u0 = besselj(0,n*x);
% end
%
%% Code Boundary Conditions
% Now, write a function that evaluates the boundary condition
% u\left(1,t\right)=J_0 \left(n\right)e^{-n^2 t} \ldotp
% Since this problem is in cylindrical coordinates (m = 1), pdepe automatically enforces the symmetry condition at x=0, so you do not need to specify a left boundary condition.
% 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\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 the partial derivatives of u need to be expressed in terms of the flux f\left(x,t,u,\frac{\partial u}{\partial x}\right). So the right boundary condition for this problem is
% u\left(1,t\right)-J_0 \left(n\right)e^{-n^2 t} +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] = heatbc(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] = 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
%
%% 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. Specify the points as vectors t and x. The vectors t and x play different roles in the solver. In particular, the cost and accuracy of the solution depend strongly on the length of the vector x. However, the computation is much less sensitive to the values in the vector t.
% For this problem, use a mesh with 25 equally spaced points in the spatial interval [0,1] for both x and t.

x = linspace(0,1,25);
t = linspace(0,1,25);

%
%% Solve Equation
% Finally, solve the equation using the symmetry m, the PDE equation, the initial conditions, the boundary conditions, and the meshes for x and t.

m = 1;
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

% pdepe returns the solution in a 3-D array sol, where sol(i,j,k) approximates the kth component of the solution u_k evaluated at t(i) and x(j). The size of sol is length(t)-by-length(x)-by-length(u0), since u0 specifies an initial condition for each solution component. For this problem, u has only one component, so sol is a 25-by-25 matrix, but in general you can extract the kth solution component with the command u = sol(:,:,k).
% Extract the first solution component from sol.

u = sol(:,:,1);

%
%% Plot Solution
% Create a surface plot of the solution. Since the problem is posed with cylindrical coordinates on a disc, the x values show the temperature on the disc some distance from the center, and the t values show how the temperature at a particular location changes over time.

surf(x,t,u)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])

%
% Plot the temperature change at the center (x=0) of the disc.

plot(t,sol(:,1))
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at center of disc')

%
%% 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] = 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
%----------------------------------------------

% Copyright 2019 The MathWorks, Inc.
%