# Stuck on a hideous differential equation

I'm stuck on this PDE: $$\rho \cdot C \cdot \dfrac{\partial T}{\partial \tau} + u \cdot \rho \cdot C \cdot \dfrac{\partial T}{\partial r} = \dfrac{\lambda}{r^{2}} \cdot \dfrac{\partial}{\partial r}\left(r^{2} \cdot \dfrac{\partial T}{\partial r}\right) + K - \dfrac{\rho \cdot C \cdot T}{V} \cdot \dfrac{\partial V}{\partial \tau}$$ This should describe the variation of enthalpy during a flash pyrolysis. I know all the parameters ($\rho, C, \lambda, K$) and I can obtain $u$ and $V$ from analytic expression.

Assuming I know the initial condition and two boundary conditions, I'd like to solve it numerically using matlab.

Can anybody help me out? Thanks in advance.

IC and BC: $$T_{\tau = 0} = T_{0} \hspace{2cm} 0 \leq r \leq R$$ $$\dfrac{\partial T}{\partial r}\bigg|_{r=0} = 0$$ $$\dfrac{\partial T}{\partial r}\bigg|_{r=R} = C1 \cdot (T_{m} - T_{R})$$

Yes, $\tau$ is a time variable and $r$ is the space coordinate. I'm expecting a profile of $T$ with respect to both variables.

I'm sorry if I can't provide so much useful information, but I was given a scientific paper to study where these equations (there is a lot of equations, but this one is the trickiest) are displayed but there is no mention of how they solved and got to the results.

• Can you provide a few more details, i.e. what variables you are solving for? Are you solving for $T$ over some spatial coordinate $r$? If so, what is $\tau$? And what do your boundary conditions look like? Oct 25, 2013 at 11:34
• This Question was originally posted at Math.SE, where because of its computational nature this forum was recommended. I surmise that $\tau$ is a time-like variable and $r$ is a 1D space variable. Note however the singularity of coefficients at $r=0$.
– hardmath
Oct 25, 2013 at 12:26
• Hi @user103012, and welcome to scicomp! Are you familiar with finite difference methods for PDE's?
– Paul
Oct 25, 2013 at 13:03
• You can separate this even further with $g(\tau)h(r)=T(\tau,r)$, since $\tau$ and $r$ are not coupled. Likely, the decoupled equation may have an analytical solution.
– Bort
Oct 28, 2013 at 0:31

It's just an advection diffusion equation in one space dimension. I don't know how to do it in matlab, but there are tutorial program in most of the open source finite element libraries (deal.II, fenics, libmesh, ...) that demonstrate how to solve such problems. Solving this equation with a finite element code isn't going to be particularly difficult as long as either (i) the diffusion dominates the advection, or (ii) you properly stabilize the equation, or (iii) you can afford a sufficiently fine mesh, which you likely can because your model is only 1d.

(Disclaimer: I'm one of the authors of deal.II.)

• I'm pretty sure that your library would solve this and I'd gladly give it a try but I don't think I have the competence to write down this equation in a C++ code. Oct 25, 2013 at 12:51
• There is a singularity in $r=0$ in the coefficients of the equations, so it's not exactly the classical AD equation. Oct 25, 2013 at 13:17
• @Dr_Sam: The BC $\partial_rT=0|_{r=0}$ should prevent that singularity, no? Oct 26, 2013 at 3:15
• The singularity disappears if you take the appropriate inner product for spherically symmetric problems, which is $(u,v)=\int_0^R u(r) v(r) r^2 \; dr$. Oct 26, 2013 at 10:40

Given the way they the boundary conditions are given I would suggest a finite difference code as hinted by Paul up above. It seems that you can set up a grid for $\tau$ and $r$ and pick a forward centered scheme for $\tau$ and possibly a centered difference scheme for $r$. I can go into more detail if you know finite differencing, of if you eventually want to use this method. Just as a first step though you may want to get all of the derivative with respect to to $r$ on one side and those w.r.t to $\tau$ on the other.

I know that Matlab has built in ODE solvers but I'm not sure about PDE solvers. I think you will learn the most by making a first attempt by yourself, by writing some code in a .m file.

• Thanks @Paul. You see, I'm not familiar with PDE or FEA. I searched for some tutorial of how to solve PDEs in matlab but I've never been able to find an example that suits my equation. Oct 26, 2013 at 8:51

I'd also suggest a finite difference code (as suggested by Paul and AshikIdrisy), but rather than discretizing in both $r$ and $\tau$, I would only discretize in $r$. Semi-discretizing the differential equations in $r$ will yield a system of ODEs in $\tau$, which can then be solved in MATLAB using one of MATLAB's built-in ODE solvers.

As far as PDE solvers go, you could try using MATLAB's pdepe solver, since your equation looks parabolic.