Sample Problems > Usage > Misc > data

{ DATA_FITTING.PDE

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'

coordinates cartesian2

variables p(Threshold = 0.001) ! the raw data field

global variables

A(Threshold = 0.001) ! The Fit amplitude

W(Threshold=0.001) ! The Fit half-width

definitions

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

initial equations

A : integral(2*fitg*dfdk) = 0

W: integral(2*fitg*dfdw) = 0

equations

p : dt(p) =div(k*grad(p))

A : integral(2*fitg*dfdk) =0

W: integral(2*fitg*dfdw) =0

boundaries

region 1

start(-box,-box)

line to (box,-box) to (box,box) to (-box,box) to close

initial values

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

plots

for cycle=1

elevation(fitf,p) from (-box,0) to (box,0)

report(A) report(integral(2*fitg*dfdk))

contour(fitf)

contour(p)

contour(fitg) as "Fit Error"

contour(dfdk)

contour(dfdw)

history(A)

history(W)

end

data_fitting

data_fitting