Fun With Numerical Methods in LaPlace Space

We've all solved partial diffential equations (PDE's) by the LaPlace transform method -- right? You apply the LaPlace transform to the PDE and boundary conditions by looking them up on LaPlace transform tables. After applying the transformed initial condition, you're left with an ODE that can usually be solved analytically without too much difficulty. But then comes the moment of truth -- can you invert your solution back to real space using inverse LaPlace transform tables? In a graduate-level transport class you probably won't be given a PDE that can't by inverted back with the tables, but in real-life these situations can and do arise. Plus, in real-life you usually have some data that you want to fit with your LaPlace-space model. What do you do then? Well, if you really must know...

You can make a MATLAB function that numerically inverts your model and use the function nlinfit to do a regression of your function to the data! Sure, you could do this on FORTRAN with IMSL routines, but 'in this economy' we don't think in terms of punchcards. Plus, Matlab has built in non-linear regression functions. The first thing you will need is a good numerical LaPlace inversion function for MATLAB -- we use invlap.m, which can be downloaded here. In fact, we recommend going through this extensive tutorial on how invlap.m is employed before attempting to write your own code. Basically, you write one function that accepts parameter values and returns the LaPlace space solution only as a function of s, the LaPlace space variable. Then you write another function in the form of "ordinate = model(parameters,abcissa)", in which parameters will be the parameters you want fitted, and abcissa will be the x-values of your data. This model should be able to take an arbitrary vector of parameters and vector of x-values and return a vector of y-values by calling invlap. Try to guess and check your parameters vector so it matches up with your data visually well. Then, you can start the fitting with something like [parameters,resid,J,Sigma] = nlinfit(abcissa,ordinate,'model',parameters) and even get error bars on the fitted parameters with nlparci.

Of course, things may not be as simple as we've briefly described. You can't employ constraints with nlinfit, so your initial guesses have to be pretty good in order for the regression to converge. When our data sets had 300+ data points, the fitting process was prohibitively slow, so we actually had to arbitrarily compress our data by deleting superfluous points. Also, for our five-parameter model, one of the parameters was consistently not optimized by nlinfit when all five-parameters were attempted to be fit simultaneously. So, we had to iteratively alternate between doing two-parameter and five-parameter fits in order to truly optimize all of the points. Our next task should be to write a script that does this automatically. But, the reward has been being able to apply a LaPlace space model to real-space data in an optimized, satisfying way. We'll be sure to update this post with a link to the code when it is published.