plot and fit surface

Hello everyone,

I hope that someone can help me.

I'm trying to plot several points in the space and to fit them so I can get an mathematical formula like this topic https://de.mathworks.com/help/curvefit/polynomial.html#bt9ykh7 "Fit and Plot a Polynomial Surface"

I want to get a quadratic polynomial so I'm trying to use 'poly22' function. This is my code:

X = [0; 1; 1; -1; -1; 0.3] 
Y = [0; 1 ; -1; -1; 1; 0.5]
Z = [0.9; 0.3; 0; 0.2; 0.6; 1]
plot3(X,Y,Z,'or')
z = Z;
fitsurface=fit([X,Y],z, 'poly22','Normalize','on')
plot(fitsurface, [X,Y],z)

I obtained this results:

Linear model Poly22:

fitsurface(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2

where x is normalized by mean 0.05 and std 0.9028

and where y is normalized by mean 0.08333 and std 0.9174

Coefficients:

p00 = 0.9097

p10 = -0.2332

p01 = 0.2645

p20 = -1.07

p11 = -0.02071

p02 = 0.5786

But I am not sure about them and the syntax.

Could you help me please?

Thank you a lot!

Laura

Answer 1

[Mahesh Taparia]

Hi

The syntax and code which you have used is correct. The above code will fit a 2nd degree polynomial to the given data points and this is what you need.

[laura bagnale]

Thank you very much Mahesh for your kind help and reply.

I would like to ask you another help if you are available.

If I wanted to find a surface in space passing through 12 points and then find the interpolating polynomial function of second degree in x y and z (w = f(x,y,z)), should I use poly222? Is the procedure the same or do you have some suggestions?

Thank you very much!

Laura

[Walter Roberson]

There is no poly222 model.

fit() is not designed for three independent variables and one dependent variable.

You would probably be better off forming a vandermode matrix and using \ to do the fitting.

let's see, do you have enough points?

p000 + p001*z + p010*y + p100*x + p002*z.^2 + p020*y.^2 + p200*x.^2 + p011*y.*z + p101*x.*z + p110*x.*y

10 coefficients, assuming that "second degree" does not have x*y*z or a square times a different variable (total degree for term no more than 2)

10 coefficients and 12 samples is mathematically possible. Just do not expect the interpolation to be mathematically stable: a bit of noise could affect the fit a fair bit.

[laura bagnale]

Thank you very much, Walter!

This is very helpful!

I will do as you suggest and hope to be able to do so! :)

Is this the only matlab resource about this topic? https://it.mathworks.com/help/matlab/ref/vander.html

Thank you very much again!

Laura

[Walter Roberson]

You will probably need to construct the matrix by hand. Assuming column vectors and r2015b or newer,

A = [ones(size(x)), x.^(1:2), y.^(1:2), z.^(1:2), x.*y, x.*z,  y.*z]
b = w;

p = A\b;

order of results:

p000, p100, p200, p010, p020, p001, p002, p110, p101, p011

You can rearrange the columns if you prefer.

[laura bagnale]

Thank you very much fo your answer!

I did the following:

x = [0; 1; 1; -1; -1; 0.3; 0.6; 0.2; 0.4; 0.3; 0.1; 0.5]
y = [0; 1; -1; -1; 1; 0.5; 0.2; 0.4; 0.7; 0.8; 0.9; 0]
z = [0.9; 0.4; 0; 0.1; 0.3; 0.8; 0.5; 0.4; 0.3; 0.7; 0.2; 0.1]

By assuming these 12 points (let's say experimental points where x, y and z are independent variables) I constructed the matrix as you suggested

A = [ones(size(x)), x.^(1:2), y.^(1:2), z.^(1:2), x.*y, x.*z, y.*z]

That represents the Vandermonde matrix, right?

I'm sorry but I didn't understand what do you mean with b = w

b = w;
p = A\b;

How can I assume b = f(x,y,z)=w?

Then I would like to fit them by a hypersurface like f(x,y,z).

Could you still help me, please?

Thank you very much!

Laura

[Walter Roberson]

In order to be able to fit w = f(x,y,z) by a multinomial, then you have to have known w values.

After having constructed the vandermode ( A ) then you would have a classic linear system, that A*p=w which would be solved by p = A\w

But that depends upon having known w.

[laura bagnale]

Thank you very much again for your support!

I have done this with known w-values and it seems to work well.

This is the code

x = [0; 1; 1; -1; -1; 0.3; 0.6; 0.2; 0.4; 0.3; 0.1; 0.5]
y = [0; 1; -1; -1; 1; 0.5; 0.2; 0.4; 0.7; 0.8; 0.9; 0]
z = [0.9; 0.4; 0; 0.1; 0.3; 0.8; 0.5; 0.4; 0.3; 0.7; 0.2; 0.1]
w = [0.3; 1; 0; -1; 0.5; 0.7; 0.5; 0.6; 0.9; 0.2; 0.3; 0.1]


A = [ones(size(x)), x.^(1:2), y.^(1:2), z.^(1:2), x.*y, x.*z,  y.*z]
b = w;
p = A\b;

Now I have the A matrix and the p-coefficients vector.

Once I have the function in three variables, how can I fit it?

I used to use fitsurface for functions of two variables, but now I have the function w=f(x,y,z) how can I represent it graphically?

cftool? fimplicit3?

Thank you

[Walter Roberson]

format long g
x = [0; 1; 1; -1; -1; 0.3; 0.6; 0.2; 0.4; 0.3; 0.1; 0.5];
y = [0; 1; -1; -1; 1; 0.5; 0.2; 0.4; 0.7; 0.8; 0.9; 0];
z = [0.9; 0.4; 0; 0.1; 0.3; 0.8; 0.5; 0.4; 0.3; 0.7; 0.2; 0.1];
w = [0.3; 1; 0; -1; 0.5; 0.7; 0.5; 0.6; 0.9; 0.2; 0.3; 0.1];


A = [ones(size(x)), x.^(1:2), y.^(1:2), z.^(1:2), x.*y, x.*z,  y.*z]

A = 12×10

1 0 0 0 0 0.9 0.81 0 0 0 1 1 1 1 1 0.4 0.16 1 0.4 0.4 1 1 1 -1 1 0 0 -1 0 0 1 -1 1 -1 1 0.1 0.01 1 -0.1 -0.1 1 -1 1 1 1 0.3 0.09 -1 -0.3 0.3 1 0.3 0.09 0.5 0.25 0.8 0.64 0.15 0.24 0.4 1 0.6 0.36 0.2 0.04 0.5 0.25 0.12 0.3 0.1 1 0.2 0.04 0.4 0.16 0.4 0.16 0.08 0.08 0.16 1 0.4 0.16 0.7 0.49 0.3 0.09 0.28 0.12 0.21 1 0.3 0.09 0.8 0.64 0.7 0.49 0.24 0.21 0.56

b = w;
p = A\b;
pcell = num2cell(p);
[p000, p100, p200, p010, p020, p001, p002, p110, p101, p011] = pcell{:};
f = @(x,y,z) p000 + p001*z + p010*y + p100*x + p002*z.^2 + p020*y.^2 + p200*x.^2 + p011*y.*z + p101*x.*z + p110*x.*y

f = function_handle with value:
    @(x,y,z)p000+p001*z+p010*y+p100*x+p002*z.^2+p020*y.^2+p200*x.^2+p011*y.*z+p101*x.*z+p110*x.*y

and now f is the fitted function. You do not need to do any more fitting on it.

Representing it graphically is more of a challenge.

N = 20;
xvec = linspace(min(x), max(x), N);
yvec = linspace(min(y), max(y), N+1);
zvec = linspace(min(z), max(z), N+2);
[X,Y,Z] = meshgrid(xvec, yvec, zvec);
W = f(X, Y, Z);
L = [-.7, -.5, -.3, 0, .3, .5, .7];
for K = 1 : length(L)
  isosurface(X, Y, Z, W, L(K));
end
legend(string(L))

[laura bagnale]

I would really like to thank you Walter, for helping me with my question.

Now everything is clearer and I have learned important things from your answers, especially from the last code.

I will try to do it again to fix all the situation!

Thank you so much for your support!

Laura

[Walter Roberson]

By the way, the reason I used N, N+1, N+2 for the sizes, is that it is easy to get the order of parameters and their orientation confused for some of the graphics operations such as surf(), so I have the habit of deliberately making the sizes different so that when I look at size() of the arrays and they only match up one way.

[laura bagnale]

Thank you very much for clarifying this further detail to me!

[laura bagnale]

Hi,

if you have time, could you help me a little bit again, please?

I wanted to minimize the previous function and to see the minimum point on the graph. Is my following code correct?

fun = (@(x,y,z) -0.3973 + 2.8392*z + 0.5153*y + 0.7036*x  -2.2113*z.^2 -0.1385*y.^2 + 0.3061*x.^2 - 0.5117*y.*z - 1.3679*x.*z - 0.0074*x.*y);

p = optimproblem
p.Description = "minimization"
x = optimvar("x", "LowerBound", -1, "UpperBound",1);
y = optimvar("y", "LowerBound", -1, "UpperBound",1);
z = optimvar("z","LowerBound", -1, "UpperBound",1 );

p.Objective = -0.3973 + 2.8392*z + 0.5153*y + 0.7036*x  -2.2113*z.^2 -0.1385*y.^2 + 0.3061*x.^2 - 0.5117*y.*z - 1.3679*x.*z - 0.0074*x.*y
initialPt.x = 0.1;
initialPt.y = 0.1;
initialPt.z = 0.1;
hold on

[sol,fval,exitflag,outout] = solve(p,initialPt)
plot3(sol.x,sol.y, sol.z, 'or','LineWidth',4)
hold off

Is there a more correct way to do so? I tried with these cases:

https://it.mathworks.com/help/matlab/ref/fminsearch.html

https://it.mathworks.com/matlabcentral/answers/517666-how-to-find-the-minimum-of-a-three-variable-function-optimization

But without success.

Thank you very much!

Laura

[Walter Roberson]

fun = @(x,y,z) -0.3973 + 2.8392*z + 0.5153*y + 0.7036*x  -2.2113*z.^2 -0.1385*y.^2 + 0.3061*x.^2 - 0.5117*y.*z - 1.3679*x.*z - 0.0074*x.*y;

syms x y z
w = fun(x,y,z)

w = 

solx = solve(diff(w,x),x)

solx = 

w2 = subs(w, x, solx)

w2 = 

soly = solve(diff(w2,y),y)

soly = 

w3 = subs(w2, y, soly)

w3 = 

solz = solve(diff(w3,z), z)

solz = 

Z = solz

Z = 

Y = subs(soly, z, Z)

Y = 

X = subs(subs(solx, y, soly), z, Z)

X = 

Xd = double(X)

Xd = 0.0396

Yd = double(Y)

Yd = 0.8851

Zd = double(Z)

Zd = 0.5273

and your task now is to determine whether (X, Y, Z) is the maxima or minima .

[laura bagnale]

Ok, thank you a lot!

It's a bit complicated for me I guess, but I'll study your code and I'll try to do my best!

Laura

[laura bagnale]

I have to conclude that my code was wrong for my purpose, right?

Thank you!

[Walter Roberson]

When you have a minimization problem, analytic solutions are best unless they would take an undue amount of time.

(There are some minimization problems that can be approached probabilisticly to get a likely solution in a relatively short time, but proving that the answer is the best possible might take a long time. There is a famous mathematical problem involving one of the largest numbers ever invented, literally too large to write down in this universe... for a situation where it is suspected that the real minimum is 6. So sometimes it really does not pay to do a complete analysis. But in a situation like the function you have, you might as well go for the analysis and so be sure that you have the right solution.

[laura bagnale]

I understand.

Thank you very much for the quick answers, the explanation and the support!

Kind Regards,

Laura