Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Suggestions for scheme to use to solve PDE numerically

  1. Feb 20, 2012 #1
    Hello everyone,
    I am trying to model the process of laser ablation on a material using MATLAB. The governing equation is of the form:
    ∂T(x,t)/∂t = ∂/∂x(A*∂T/∂x) + B*exp(-C*t2)*exp(-D*x)
    with one Initial condition and two boundary conditions.
    Using the built-in 'pdepe' function in Matlab gave inaccurate results, so I have been reading material on solving it numerically. This being a non-linear parabolic pde (correct me if wrong), I am facing difficulty on how to proceed. All the materials that I have referred to so far only address problems where the last term (source term) of the equation is a function of T(x,t). Please guide me.
  2. jcsd
  3. Feb 20, 2012 #2


    User Avatar
    Gold Member

    Is A a function of x or t or T?

    If its not a function of any of the above this is just a non-homogenous heat equation and you can solve it analytically.

    If it's a function of x or t or both, it's still a linear equation.

    If it's a function of T then indeed this is non-linear.
  4. Feb 20, 2012 #3


    User Avatar
    Gold Member

    Double Post. Sorry.
  5. Feb 20, 2012 #4
    Usually for these problems the Crank-Nicolson method suffices. Note that your source term is quite sensitive to t and x with an exponential decay so a very fine mesh is probably needed when x<<1 and t<<1 to capture the source term accurately.

    Also, in numerical methods, it doesn't matter (much) if the source term depends on T or not. In this case, the source term is known on every point of the grid if B,C,D are not functions of T.

    You should also be able to get the analytical solution as kai_sikorski mentioned. It is always good to have an analytical solution.
  6. Feb 21, 2012 #5
    Sorry. Forgot to mention about the coefficients earlier.
    A is a function of T, whereas B,C and D are positive constants.

  7. Feb 21, 2012 #6
    As you have noticed that i will be needing an extremely fine mesh, I was using a time step of 1ns and a 10nm step for distance x. Is it possible that using these increments might affect the stability of the solution?
    I would also like to mention that I need to study the behaviour in the first 100ns and the ablation was observed till a few μm.

  8. Feb 21, 2012 #7


    User Avatar
    Gold Member

    It might be a good idea to non-dimensionalize time and space. Quantities with units make it hard to tell what's large and what's small. Is 1ns small? Well it depends on how fast the forcing in the problem is varying. If C~(1 ns)^2, or smaller then your time step might be huge in comparison. So I would non-nimensionalize time by (C)^{-1/2}, and space by (D)^{-1}, that way the forcing is changing over order 1 quantities in space and time and the meaning of what a small time step is will be more clear.

    Stability is usually related to the time step being too large compared to the space step. I just read the documentation for pdepe though and it seems that it actually chooses an appropriate time step on it's own, the time points you specify are just the time points it will give you a solution at, but they have nothing to do with what time step it actually used. So presumably it chooses something that should be stable.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook