MATLAB Answers

How to fit data to the partial differential equation for 1-D diffusion?

5 views (last 30 days)
Archishman Ghosh
Archishman Ghosh on 13 May 2021
Edited: Archishman Ghosh on 14 May 2021
I have some data represented as a surface below with time and distance axes. I would like to fit this to all possible solutions for the diffusion equation that is:
where S is the signal representing concentration here, t is time and x is distance.
I think the pdepe can do this but I'm not sure how to implement it? I'm familiar with fitting data to a function which often is the analytical solution to diff. eqns, but am not sure how to solve and fit at the same time. Could someone explain/demonstrate this to me?
Example matrix data attached. All data(:,1) are on the distance axis and data(1,:) are on the time axis.
Thanks!
  4 Comments
Archishman Ghosh
Archishman Ghosh on 14 May 2021
Yes. There are 84 traces (each one a frame in time) that are 45 distance units long each. Each trace shows how the signal(=concentration) is distributed over the 0 to 44 distance units. As the values at the edges are different over each frame going from 0 to 83, the boundary conditions also evolve. Thank you so much for your interest.

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 14 May 2021
Edited: John D'Errico on 14 May 2021
Ok. I had some time to look at your surface more closely, It does not fit that PDE. It cannot do so. Now that I have your data...
S is a 45x84 array. Given your plot label,s, time moves from 0-83, and distance from 0-44. (Pick your own units for those variables.)
size(S)
ans =
45 84
plot(0:44,S(:,1))
xlabel 'Distance'
I''ve plotted S, as a function of distance at time 0. You want to use a diffusion model on this system, but have posed no forcing term in the equation of any form. Therefore, since this is a diffusion process, as time moves foward, this system will approach an AVERAGE level. Diffusion is an averaging thing. The level at a long time in the future, IF that model makes sense at all,
mean(S(:,1))
ans =
1559.1
would be somethign near 1559.
But it CANNOT move to a much higher average level, unless your model includes a forcing term of some sort, where there is addtional "stuff" being injected over time. There is no such term in the model you posed.
hold on
plot(0:44,S(:,end),'r-')
So the line in red is the process after a relatively long time, the blue line is the process at time 0. Sorry, but this system is NOT a simple diffusion process. You CANNOT fit that model to this surface. The two are completely inconsistent.
Unless you have a better model for this process, trying to do the fit you asked to do would be meaningless here...
===================================
Edit:
I ran out of time last night to finish what I wanted to say. So now that we know the PDE that you posed CANNOT possibly solve the problem, is there a simple variation thereof that would do so? This now becomes entirely a problem in mathematical modeling, not for MATLAB, at least not immediately. Once we know what mathematics we need to use, then we can come back to MATLAB to write code.
Looking at the surface plot you show, we see the spatial variation is actually pretty small. Far more so than the temporal variation. So, temporarily, I'll assume that this is purely a temporal problem, so a 1-d ODE. What ODE has a solution that looks like an exponential rise to a constant asymptote? That one is easy.
dS/dt = -k*S + u
Here we don't know the rate constant. It might even be a function of x. But if we integrate that ODE, we would get a solution that behaves like your process.
syms S(t) k u
dsolve(diff(S,t) == -k*S + u)
ans =
(u - C1*exp(-k*t))/k
Here we would see a process that rises to an upper asymptote of u. The rate constant for the exponential is k. And C1 is an undetermined constant.
So now we can put this back into a PDE setting. It looks like after a long time, the process rises to an upper asymptote, one that is the same for all x. But the rate constant may vary. So I'll assume that k is an unknown function of x, though for a first pass, I might assume it is truly constant. Anyway, the PDE now becomes...
St = -k(x)*S + u + D*Sxx
Here St is the first partial of S with respect to t, and Sxx is the second partial with respect to x.
If k(x) is truly a constant, then the PDE becomes
St = -k*S + u + D*Sxx
This is still NOT something we can throw into PDEPE. But, since we have data, we can now use that data.
We can approximate those partial derivatives using finite differences. Then with k, u, and D as unknowns, the problem will become a simple linear least squares. In fact, we could then use backslash to estimate the unknown constants.
The next question is, what finite difference approximations would we employ? (By the way, all of this is enough that I could probably write a paper on the modeling, analysis and solution.) But your data is noisy. Enough so that a simple finite difference will be far more noisy. So I might use a trick from Savitsky-Golay now. For example, perhaps we might find a 7 point finite difference approximation to a series for the first derivative at any point, assuming the function is locally quadratic.
X = (-3:3)';
A = [ones(7,1),X,X.^2];
Ap = pinv(A);
StFDA = Ap(2,:)
StFDA =
-0.10714 -0.071429 -0.035714 2.1371e-17 0.035714 0.071429 0.10714
The vector StFDA represents the coefficients of the best finite difference estimator of the first derivative of a vector, subject to noise, and assuming the process is locally quadratic. We can apply that finite difference approximation to the surface array using conv2.
Similarly, we can find another such approximation for the second partial derivative of the surface, though I'll need a factor of 2 to get the necessary coefficients.
X = (-3:3)';
A = [ones(7,1),X,X.^2];
Ap = pinv(A);
SxxFDA = Ap(3,:)*2
SxxFDA =
0.11905 2.4547e-17 -0.071429 -0.095238 -0.071429 -3.2735e-18 0.11905
I know, you still don't see where this is going. But I do, and since this is an interesting problem in modeling, press forwards...
St_approx = conv2(S,StFDA','same');
Sxx_approx = conv2(S,SxxFDA,'same');
% trim them, because these approximations are not valid around the edges
St_approx_trim = St_approx(3:end-2,3:end-2);
Sxx_approx_trim = Sxx_approx(3:end-2,3:end-2);
S_trim = S(3:end-2,3:end-2);
kuD = [-S_trim(:),ones(numel(S_trim),1),Sxx_approx_trim(:)]\-St_approx_trim(:)
kuD =
0.11611
513.03
0.7523
So my best estimate of the rate constant k is 0.11611. The estimated upper asymptote is u/k,
kuD(2)/kuD(1)
ans =
4418.4
and my estimate for D is 0.7523. Note that the upper asymptote is not too bad, compared to the data.
plot(0:44,S(:,1))
hold on
plot(0:44,S(:,end),'r-')
yline(4418.4);
We see the upper asymptote does seem to pretty much pass through that final time point.
With a little more effort, I could have chosen a better form for k(x), if we wanted it to be non-constant. (k really does seem to vary with x, if I look at that surface.) But what I've already done is probably way more work than you ever thought was going to be necessary.
Finally, now that we know the estimates of k, u, and D, if we choose to provide intelligent boundary conditions and boundaty values on the solution, we could even throw this into PDEPE. The resulting prediction would now pretty nicely approximate the surface plot from your data.
  1 Comment
Archishman Ghosh
Archishman Ghosh on 14 May 2021
Thank you for your prompt reply. You are absolutely correct to say that it will not satisfy the differential equation with a constant total concentration. I should have added, particles whose concentration is represented by S are indeed being injected (or rather are diffusing in) from the outer boundary provided here. Essentially S is the concentration of particles inside a circular compartment wherein the same particles are diffusing across a boundary to replenish it. The dimension has been reduced to 1-D by radially averaging the S values inside that spherical compartment.
So my understanding is, one has to update the boundary conditions for every single single trace from 0 to 83 and solve the equation for each trace. Indeed you plotted the book ends to the traces in your image and I've attached the complete set.
And the boundary conditions would evolve as shown in the second image.
I'm also attaching a cartoon describing my problem. I'm better with pictures than words.
Edit: I'm also attaching representative images of the real data from which the traces are plotted. These are initial experiments done on a slower instrument. Getting smoother data won't be an issue. The analysis part is more challenging to me.
Time = 0
Time = 4
Time = 30

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!