Hello,
I’m a PhD student in modeling mechanical behavior of materials. For my work, i need to fit a surface on a 3d curve x,y,z. I had a look in all libraries but i had difficulties to find an understandable working example. The best code i have found is the one from Alan Miller : lsp and quadsurf but the resulting surface does not fit my input data for some reason…
Does anyone know where i can find a code to do this 3d fitting with polynomials?
Thank you for your help.
Romain
Welcome to the forum. Do you have a working code for multiple regression (linear least squares)? (I have used Burkardt’s qr_solve.) If so, using it to fit a quadratic surface in two dimensions is simple – something like this:
program xsurface
integer, parameter :: nobs = 1000, dp = kind(1.0d0), ncoeff = 6
real(kind=dp) :: x(nobs),y(nobs),z(nobs),a(nobs,ncoeff),coeff(ncoeff)
! read x, y, z
! construct design matrix A
a(:,1) = 1.0_dp ! intercept
a(:,2) = x
a(:,3) = y
a(:,4) = x**2
a(:,5) = y**2
a(:,6) = x*y
call fit(a,z,coeff) ! need a subroutine that solves A*coeff = Z in the least-squares sense
end program xsurface
I read this to mean that you have one surface to fit.
I like to use gnuplot for this. The fit() function uses the nonlinear least-squares Marquardt-Levenberg algorithm. It is described in the documentation and there are examples at gnuplot demo script: fit.dem It isn’t Fortran but it does have good graphics, which are important to assess your result.
In simple cases you can use the LINEST function in Excel to do a multiple linear regression using the values of the polynomials as the independent variables.
Thank you Beliavsky, the code you suggest adapted a bit with the qr_solve is working !! 
You made my day !
Thank you
The LM method is for nonlinear least-squares. The problem posed by the OP is a linear least-squares problem.
Sorry, that does not make sense. You can find an infinite number of surfaces on which a given curve can be located. Consider, for example, a straight line in 3-D space. If you find one plane that contains that line, you can rotate the plane about the line and all the rotated positions also “fit” the line.
Yes. LM is overkill here, but it works well provided it is given a reasonable starting point. I like it because it is allows to nonlinear approximations.
Or, just use LAPACK then? Seems to be the straight-forward way.
Yes. The Lapack95 subroutine call is simply
call gels(a,z)
if a suitable USE statement with the interface for GELS has been provided in the caller. The solution vector is returned in z.