This example uses GLOBAL VARIABLES to form a least-squares fit to a Gaussian data distribution,
and then follows the fit as the Gaussian diffuses out.
The basic process of least-squares fitting seeks to minimize the integral of the square of the
difference between the given data and the analytic fit function:
minimize G = Integral (( F - P)^2 * dV)
where F is the analytic fit and P is the array of given data.
The technique is to find a stationary point in the derivatives of G with respect to the fit parameters.
In our case, we choose F(x,y) = A*exp(-R^2/W^2), where A is the amplitude and W is the half-width
of the fitted Gaussian, and R is the radius sqrt(x^2+y^2), (a pre-defined name in FlexPDE).
With this definition, we can define the fit equations
dG/dA = Integral(2*(F-P)*dF/dA) = 0
dG/dW = Integral(2*(F-P)*dF/dW) = 0
We start by solving the fit equations in an INITIAL EQUATIONS section, then proceed to solve the
fit equations simultaneously with the diffusion of the data field.
title 'Least-Squares Data Fitting'
variables p(Threshold = 0.001) ! the raw data field
A(Threshold = 0.001) ! The Fit amplitude
W(Threshold=0.001) ! The Fit half-width
box = 1.5 ! The Domain Size
wdata = 0.1 ! The half-width of the data Gaussian
k = 0.1 ! the diffusivity of the data equation
! Force denser mesh at the peak of the data
mesh_den = 50
mesh_density = mesh_den*exp(-(x^2+y^2)/wdata^2)
! The Gaussian fit equation terms
fitf = A*exp(-r^2/W^2)
dfdk = fitf/A
dfdw = fitf * (2*r^2/W^3)
fitg = fitf - p
A : integral(2*fitg*dfdk) = 0
W: integral(2*fitg*dfdw) = 0
p : dt(p) =div(k*grad(p))
A : integral(2*fitg*dfdk) =0
W: integral(2*fitg*dfdw) =0
line to (box,-box) to (box,box) to (-box,box) to close
p = exp(-(x^2 + y^2)/wdata^2)
A = 0.9 ! slightly erroneous first guess amplitude
W= 1.1*wdata ! slightly erroneous first guess half-width
time 0 to 1
elevation(fitf,p) from (-box,0) to (box,0)
contour(fitg) as "Fit Error"