Efficient way to apply fsolve to large 2D array element

Hello

I have 2 nonlinear equations. I figure out that it can be solved by fsolve. The outputs of the 2 nonlinear equations are the [x,y] position of a point. The inputs of the system are initial estimate of the point [x0,y0] and some constants. I first tried to solve the nonlinear equations for one point, with the following code:

options = optimoptions('fsolve', 'Display', 'off');
x0 = [ud(1,1);vd(1,1)]; % Initial guess of undistorted point [x,y]
[x, fval] = fsolve(@myfunction,x0,options, kc, x0);
function F = myfunction(x, kc, x0)
r = sqrt(x(1)^2 + x(2)^2);
F = [ -x0(1) + (1 + kc(1)* r.^2 + kc(2)* r.^4 + kc(5)* r.^6 ).*x(1) + 2*kc(3).*x(1).*x(2) + kc(4)*(r.^2+2*x(1).^2), ...
      -x0(2) + (1 + kc(1)* r.^2 + kc(2)* r.^4 + kc(5)* r.^6 ).*x(2) + kc(3)*(r.^2 + 2*x(2).^2) + 2*kc(4) .*x(1) .*x(2)];
end

The above code runs perfectly fine with one point [x,y]. However, I would like to apply it to large set of points. The set of points are stored in two m x n arrays, ud and vd, where [ ud(row,col), vd(row,col) ] forms the initial estimiate [x0,y0] for each array element. I tried to use nested for loop to run the above fsolve code for each array element. While it works, it takes decades to complete because m and n are large (>1500x1500). Is there any efficient way to apply the above code to each array element? Or is there any other way to do the same thing? Thanks!

[Bruno Luong]

Sorry for a stupid question: but you have 2 unknowns and 2 equations. They look reasonable, so the number of solutions should be small (no oscillating function involve). Why you look at solutions with 1500^2 stating points? What is the goal? They must converge to few solutions of the equation IMO.

[Yuk Ching Leung]

Hi For each element in the 2D array, there are 2 nonlinear equations and 2 unknowns. And those equations and unknowns are independent of each other. Therefore, there are essentially >1500 set of 2 nonlinear equations to solve.

[Bruno Luong]

I see, so x0 is not only "forms the initial estimate" but also parametrize the function F.

[Walter Roberson]

Dang, I had a solution written up but accidentally Cancel instead of Submit :(

Summary: there is a solution

x2 = roots() of a degree 9 polynomial

after which you select the real-valued x2, then

x1 = (4.*ks3.*ks4.*ks5.*(ks3.^2+ks4.^2).^2.*x2.^7+ks5.*((ks1.*x01+ks4).*ks3.^4+2.*x02.*x01.*ks2.*ks3.^3+((-x01.^3+3.*x01.*x02.^2).*ks5+2.*ks4.^3+2.*ks1.*ks4.^2.*x01+(ks2.*x01.^2-ks2.*x02.^2).*ks4).*ks3.^2+2.*((3.*x01.^2-x02.^2).*ks5+ks2.*ks4.*x01).*x02.*ks4.*ks3+ks4.^2.*((x01.^3-3.*x01.*x02.^2).*ks5+ks4.*(ks1.*ks4.*x01+ks2.*x01.^2-ks2.*x02.^2+ks4.^2))).*x2.^6+((4.*ks2.*ks4-ks5.*x01).*ks3.^5+3.*ks3.^4.*ks4.*ks5.*x02+(4.*ks2.*ks4.^3-10.*ks4.^2.*ks5.*x01).*ks3.^3+10.*ks3.^2.*ks4.^3.*ks5.*x02-ks3.*ks4.^4.*ks5.*x01-ks4.^5.*ks5.*x02).*x2.^5+((-ks5.*x01+ks2.*(ks1.*x01+ks4)).*ks3.^4+2.*x02.*(ks2.^2.*x01+ks4.*ks5).*ks3.^3+((-3.*ks4.^2.*x01+(-3.*ks1.*x01.^2+ks1.*x02.^2).*ks4-ks2.*x01.^3+4.*ks2.*x01.*x02.^2).*ks5+ks2.*ks4.*(ks1.*ks4.*x01+ks2.*x01.^2-ks2.*x02.^2+ks4.^2)).*ks3.^2+4.*((-(1./2).*x01.^3+(1./2).*x01.*x02.^2).*ks5+ks4.*(ks1.*ks4.*x01+(1./2).*ks2.*x01.^2-(1./2).*ks2.*x02.^2+ks4.^2)).*x02.*ks5.*ks3+((3.*x01.^2-x02.^2).*ks5+ks2.*ks4.*x01).*x02.^2.*ks5.*ks4).*x2.^4+4.*((ks1.*ks4-(1./4).*ks2.*x01).*ks3.^4+(3./4).*x02.*(ks2.*ks4-ks5.*x01).*ks3.^3-(1./4).*((x01.^2-5.*x02.^2).*ks5+ks2.*ks4.*x01).*ks4.*ks3.^2-(1./4).*ks4.^2.*x02.*(ks2.*ks4-3.*ks5.*x01).*ks3-ks4.^3.*ks5.*x02.^2).*ks3.*x2.^3+(8.*ks3.^6.*ks4+(ks1.*ks4+x01.*(ks1.^2-ks2)).*ks3.^4+2.*x02.*(-ks5.*x01+ks2.*(ks1.*x01+ks4)).*ks3.^3+2.*((ks1.*x01+(3./2).*ks4).*ks5+(1./2).*ks2.^2.*x01).*x02.^2.*ks3.^2+2.*ks2.*ks5.*x01.*x02.^3.*ks3+ks5.^2.*x01.*x02.^4).*x2.^2+2.*((ks1.*x01+(3./2).*ks4).*ks3.^3-(1./2).*x02.*(ks1.*ks4-3.*ks2.*x01).*ks3.^2-x02.^2.*(ks2.*ks4-2.*ks5.*x01).*ks3-(3./2).*ks4.*ks5.*x02.^3).*ks3.^2.*x2+ks3.^6.*x01-3.*ks3.^5.*ks4.*x02)./(2.*ks5.*(ks3-ks4).*(ks3+ks4).*(ks3.^2+ks4.^2).^2.*x2.^6+ks5.*(ks3.^5+ks1.*ks3.^4.*x02+(-ks2.*x01.^2+ks2.*x02.^2+2.*ks4.^2).*ks3.^3+2.*((-(3./2).*x01.^2+(1./2).*x02.^2).*ks5+ks1.*ks4.^2+ks2.*ks4.*x01).*x02.*ks3.^2-((2.*x01.^3-6.*x01.*x02.^2).*ks5+ks4.*(ks2.*x01.^2-ks2.*x02.^2-ks4.^2)).*ks4.*ks3+x02.*ks4.^2.*((3.*x01.^2-x02.^2).*ks5+ks1.*ks4.^2+2.*ks2.*ks4.*x01)).*x2.^5+2.*(ks2.*ks3.^5+ks3.^4.*ks5.*x02+(-ks2.*ks4.^4+4.*ks4.^3.*ks5.*x01).*ks3-5.*ks4.^4.*ks5.*x02).*ks3.*x2.^4+(ks2.*ks3.^5+x02.*(ks1.*ks2+ks5).*ks3.^4+((-2.*ks4.*x01+ks1.*(x01.^2+x02.^2)).*ks5-ks2.^2.*x01.^2+ks2.^2.*x02.^2+ks2.*ks4.^2).*ks3.^3+((-4.*ks1.*ks4.*x01-3.*ks2.*x01.^2+2.*ks2.*x02.^2+3.*ks4.^2).*ks5+ks2.*ks4.*(ks1.*ks4+2.*ks2.*x01)).*x02.*ks3.^2+4.*x02.^2.*ks5.*((-(3./4).*x01.^2+(1./4).*x02.^2).*ks5+ks1.*ks4.^2+ks2.*ks4.*x01).*ks3+ks4.*ks5.*x02.^3.*(ks2.*ks4+4.*ks5.*x01)).*x2.^3+2.*(ks1.*ks3.^4+ks2.*x02.*ks3.^3+(((1./2).*x01.^2+x02.^2).*ks5-ks1.*ks4.^2+2.*ks2.*ks4.*x01).*ks3.^2-3.*ks4.*x02.*(ks2.*ks4-ks5.*x01).*ks3-(9./2).*ks4.^2.*ks5.*x02.^2).*ks3.^2.*x2.^2+(2.*ks3.^7+(-6.*ks4.^2+ks1).*ks3.^5+x02.*(ks1.^2+ks2).*ks3.^4+2.*x02.^2.*(ks1.*ks2+(1./2).*ks5).*ks3.^3+(2.*ks1.*ks5.*x02.^3+ks2.^2.*x02.^3).*ks3.^2+2.*ks2.*ks3.*ks5.*x02.^4+ks5.^2.*x02.^5).*x2+ks3.^3.*(ks1.*ks3.^2.*x02+ks2.*ks3.*x02.^2+ks5.*x02.^3+ks3.^3))

And then you look for the (x1, x2) pair that is closest to (x01, x02)

[Yuk Ching Leung]

Yes. In line with Luong's answer, this problem can be solved without applying fsolve to each array element. One thing that I would like to ask is that, what solver/function did you used to solve for the solution of x1 and x2? I tried using solve() but it doesn't seem to giving any result after a long period of runtime.

[Bruno Luong]

Just a few ideas:

What you could do is assume the inverse the same form with unknown coefs lc1,...,lc6 (instead of kci), then when combine with the direct model it gives identity map by requiring the high-order terms of the combine function = 0. This might provides a formula of lc1,...,lc6 from kc1,...,kc6 using symbolic calculation.

You might also solve numerically using FSOLVE on n>=6 on well chosen points, n not too large, then "fit" the reverse with a model to find lc1,...lc6 using the result of FSOLVE. Then use lc1,...,lc6 to compute the results for the rest of the data.

[Walter Roberson]

MATLAB unfortunately has trouble solving higher dimensional polynomials.

In my outline above,

x2 = RootOf((4*ks3^8*ks5+16*ks3^6*ks4^2*ks5+24*ks3^4*ks4^4*ks5+16*ks3^2*ks4^6*ks5+4*ks4^8*ks5)*_Z^9+(4*ks1*ks3^6*ks5*x02+8*ks1*ks3^5*ks4*ks5*x01+4*ks1*ks3^4*ks4^2*ks5*x02+16*ks1*ks3^3*ks4^3*ks5*x01-4*ks1*ks3^2*ks4^4*ks5*x02+8*ks1*ks3*ks4^5*ks5*x01-4*ks1*ks4^6*ks5*x02-4*ks2*ks3^5*ks5*x01^2+4*ks2*ks3^5*ks5*x02^2+24*ks2*ks3^4*ks4*ks5*x01*x02+8*ks2*ks3^3*ks4^2*ks5*x01^2-8*ks2*ks3^3*ks4^2*ks5*x02^2+16*ks2*ks3^2*ks4^3*ks5*x01*x02+12*ks2*ks3*ks4^4*ks5*x01^2-12*ks2*ks3*ks4^4*ks5*x02^2-8*ks2*ks4^5*ks5*x01*x02-12*ks3^4*ks5^2*x01^2*x02+4*ks3^4*ks5^2*x02^3-16*ks3^3*ks4*ks5^2*x01^3+48*ks3^3*ks4*ks5^2*x01*x02^2+72*ks3^2*ks4^2*ks5^2*x01^2*x02-24*ks3^2*ks4^2*ks5^2*x02^3+16*ks3*ks4^3*ks5^2*x01^3-48*ks3*ks4^3*ks5^2*x01*x02^2-12*ks4^4*ks5^2*x01^2*x02+4*ks4^4*ks5^2*x02^3+4*ks3^7*ks5+12*ks3^5*ks4^2*ks5+12*ks3^3*ks4^4*ks5+4*ks3*ks4^6*ks5)*_Z^8+(ks1^2*ks3^4*ks5*x01^2+ks1^2*ks3^4*ks5*x02^2+2*ks1^2*ks3^2*ks4^2*ks5*x01^2+2*ks1^2*ks3^2*ks4^2*ks5*x02^2+ks1^2*ks4^4*ks5*x01^2+ks1^2*ks4^4*ks5*x02^2+2*ks1*ks2*ks3^3*ks5*x01^2*x02+2*ks1*ks2*ks3^3*ks5*x02^3+2*ks1*ks2*ks3^2*ks4*ks5*x01^3+2*ks1*ks2*ks3^2*ks4*ks5*x01*x02^2+2*ks1*ks2*ks3*ks4^2*ks5*x01^2*x02+2*ks1*ks2*ks3*ks4^2*ks5*x02^3+2*ks1*ks2*ks4^3*ks5*x01^3+2*ks1*ks2*ks4^3*ks5*x01*x02^2-2*ks1*ks3^2*ks5^2*x01^4+2*ks1*ks3^2*ks5^2*x02^4+8*ks1*ks3*ks4*ks5^2*x01^3*x02+8*ks1*ks3*ks4*ks5^2*x01*x02^3+2*ks1*ks4^2*ks5^2*x01^4-2*ks1*ks4^2*ks5^2*x02^4+ks2^2*ks3^2*ks5*x01^4+2*ks2^2*ks3^2*ks5*x01^2*x02^2+ks2^2*ks3^2*ks5*x02^4+ks2^2*ks4^2*ks5*x01^4+2*ks2^2*ks4^2*ks5*x01^2*x02^2+ks2^2*ks4^2*ks5*x02^4+4*ks2*ks3^8+12*ks2*ks3^6*ks4^2+12*ks2*ks3^4*ks4^4+4*ks2*ks3^2*ks4^6+2*ks2*ks3*ks5^2*x01^4*x02+4*ks2*ks3*ks5^2*x01^2*x02^3+2*ks2*ks3*ks5^2*x02^5+2*ks2*ks4*ks5^2*x01^5+4*ks2*ks4*ks5^2*x01^3*x02^2+2*ks2*ks4*ks5^2*x01*x02^4-12*ks3^6*ks4*ks5*x01+8*ks3^5*ks4^2*ks5*x02-28*ks3^4*ks4^3*ks5*x01+16*ks3^3*ks4^4*ks5*x02-20*ks3^2*ks4^5*ks5*x01+8*ks3*ks4^6*ks5*x02-4*ks4^7*ks5*x01+ks5^3*x01^6+3*ks5^3*x01^4*x02^2+3*ks5^3*x01^2*x02^4+ks5^3*x02^6+2*ks1*ks3^5*ks5*x02+2*ks1*ks3^4*ks4*ks5*x01+4*ks1*ks3^3*ks4^2*ks5*x02+4*ks1*ks3^2*ks4^3*ks5*x01+2*ks1*ks3*ks4^4*ks5*x02+2*ks1*ks4^5*ks5*x01-2*ks2*ks3^4*ks5*x01^2+2*ks2*ks3^4*ks5*x02^2+8*ks2*ks3^3*ks4*ks5*x01*x02+8*ks2*ks3*ks4^3*ks5*x01*x02+2*ks2*ks4^4*ks5*x01^2-2*ks2*ks4^4*ks5*x02^2-6*ks3^3*ks5^2*x01^2*x02+2*ks3^3*ks5^2*x02^3-6*ks3^2*ks4*ks5^2*x01^3+18*ks3^2*ks4*ks5^2*x01*x02^2+18*ks3*ks4^2*ks5^2*x01^2*x02-6*ks3*ks4^2*ks5^2*x02^3+2*ks4^3*ks5^2*x01^3-6*ks4^3*ks5^2*x01*x02^2+ks3^6*ks5+3*ks3^4*ks4^2*ks5+3*ks3^2*ks4^4*ks5+ks4^6*ks5)*_Z^7+(4*ks1*ks2*ks3^6*x02+8*ks1*ks2*ks3^5*ks4*x01+8*ks1*ks2*ks3^3*ks4^3*x01-4*ks1*ks2*ks3^2*ks4^4*x02-16*ks1*ks3^4*ks4*ks5*x01*x02-32*ks1*ks3^3*ks4^2*ks5*x01^2+16*ks1*ks3^3*ks4^2*ks5*x02^2+48*ks1*ks3^2*ks4^3*ks5*x01*x02-16*ks1*ks3*ks4^4*ks5*x02^2-4*ks2^2*ks3^5*x01^2+4*ks2^2*ks3^5*x02^2+24*ks2^2*ks3^4*ks4*x01*x02+12*ks2^2*ks3^3*ks4^2*x01^2-12*ks2^2*ks3^3*ks4^2*x02^2-8*ks2^2*ks3^2*ks4^3*x01*x02-16*ks2*ks3^4*ks5*x01^2*x02+4*ks2*ks3^4*ks5*x02^3-12*ks2*ks3^3*ks4*ks5*x01^3+44*ks2*ks3^3*ks4*ks5*x01*x02^2+28*ks2*ks3^2*ks4^2*ks5*x01^2*x02-20*ks2*ks3^2*ks4^2*ks5*x02^3-4*ks2*ks3*ks4^3*ks5*x01^3-12*ks2*ks3*ks4^3*ks5*x01*x02^2+4*ks2*ks4^4*ks5*x01^2*x02+4*ks3^3*ks5^2*x01^4-12*ks3^3*ks5^2*x01^2*x02^2-32*ks3^2*ks4*ks5^2*x01^3*x02+16*ks3^2*ks4*ks5^2*x01*x02^3-8*ks3*ks4^2*ks5^2*x01^4+36*ks3*ks4^2*ks5^2*x01^2*x02^2-4*ks3*ks4^2*ks5^2*x02^4+8*ks4^3*ks5^2*x01^3*x02-8*ks4^3*ks5^2*x01*x02^3+4*ks2*ks3^7+8*ks2*ks3^5*ks4^2+4*ks2*ks3^3*ks4^4-20*ks3^5*ks4*ks5*x01+20*ks3^4*ks4^2*ks5*x02-16*ks3^3*ks4^3*ks5*x01+16*ks3^2*ks4^4*ks5*x02+4*ks3*ks4^5*ks5*x01-4*ks4^6*ks5*x02)*_Z^6+(ks1^2*ks2*ks3^4*x01^2+ks1^2*ks2*ks3^4*x02^2+ks1^2*ks2*ks3^2*ks4^2*x01^2+ks1^2*ks2*ks3^2*ks4^2*x02^2-4*ks1^2*ks3^2*ks4*ks5*x01^3-4*ks1^2*ks3^2*ks4*ks5*x01*x02^2+4*ks1^2*ks3*ks4^2*ks5*x01^2*x02+4*ks1^2*ks3*ks4^2*ks5*x02^3+2*ks1*ks2^2*ks3^3*x01^2*x02+2*ks1*ks2^2*ks3^3*x02^3+2*ks1*ks2^2*ks3^2*ks4*x01^3+2*ks1*ks2^2*ks3^2*ks4*x01*x02^2-3*ks1*ks2*ks3^2*ks5*x01^4-ks1*ks2*ks3^2*ks5*x01^2*x02^2+2*ks1*ks2*ks3^2*ks5*x02^4+4*ks1*ks2*ks3*ks4*ks5*x01^3*x02+4*ks1*ks2*ks3*ks4*ks5*x01*x02^3+ks1*ks2*ks4^2*ks5*x01^2*x02^2+ks1*ks2*ks4^2*ks5*x02^4+4*ks1*ks3^8+8*ks1*ks3^6*ks4^2+4*ks1*ks3^4*ks4^4-4*ks1*ks3*ks5^2*x01^4*x02-4*ks1*ks3*ks5^2*x01^2*x02^3+4*ks1*ks4*ks5^2*x01^3*x02^2+4*ks1*ks4*ks5^2*x01*x02^4+ks2^3*ks3^2*x01^4+2*ks2^3*ks3^2*x01^2*x02^2+ks2^3*ks3^2*x02^4+2*ks2^2*ks3*ks5*x01^4*x02+4*ks2^2*ks3*ks5*x01^2*x02^3+2*ks2^2*ks3*ks5*x02^5-4*ks2*ks3^6*ks4*x01-8*ks2*ks3^4*ks4^3*x01-4*ks2*ks3^2*ks4^5*x01+ks2*ks5^2*x01^4*x02^2+2*ks2*ks5^2*x01^2*x02^4+ks2*ks5^2*x02^6+11*ks3^6*ks5*x01^2-38*ks3^5*ks4*ks5*x01*x02-34*ks3^4*ks4^2*ks5*x01^2+27*ks3^4*ks4^2*ks5*x02^2+68*ks3^3*ks4^3*ks5*x01*x02+19*ks3^2*ks4^4*ks5*x01^2-34*ks3^2*ks4^4*ks5*x02^2-22*ks3*ks4^5*ks5*x01*x02+3*ks4^6*ks5*x02^2+2*ks1*ks2*ks3^5*x02+2*ks1*ks2*ks3^4*ks4*x01+2*ks1*ks2*ks3^3*ks4^2*x02+2*ks1*ks2*ks3^2*ks4^3*x01-ks1*ks3^4*ks5*x01^2-8*ks1*ks3^3*ks4*ks5*x01*x02-9*ks1*ks3^2*ks4^2*ks5*x01^2+9*ks1*ks3^2*ks4^2*ks5*x02^2+8*ks1*ks3*ks4^3*ks5*x01*x02+ks1*ks4^4*ks5*x02^2-2*ks2^2*ks3^4*x01^2+2*ks2^2*ks3^4*x02^2+8*ks2^2*ks3^3*ks4*x01*x02+2*ks2^2*ks3^2*ks4^2*x01^2-2*ks2^2*ks3^2*ks4^2*x02^2-10*ks2*ks3^3*ks5*x01^2*x02+2*ks2*ks3^3*ks5*x02^3-2*ks2*ks3^2*ks4*ks5*x01^3+18*ks2*ks3^2*ks4*ks5*x01*x02^2-4*ks2*ks3*ks4^2*ks5*x02^3+4*ks2*ks4^3*ks5*x01*x02^2+3*ks3^2*ks5^2*x01^4-9*ks3^2*ks5^2*x01^2*x02^2-12*ks3*ks4*ks5^2*x01^3*x02+12*ks3*ks4*ks5^2*x01*x02^3+9*ks4^2*ks5^2*x01^2*x02^2-3*ks4^2*ks5^2*x02^4+ks2*ks3^6+2*ks2*ks3^4*ks4^2+ks2*ks3^2*ks4^4-6*ks3^4*ks4*ks5*x01+6*ks3^3*ks4^2*ks5*x02-6*ks3^2*ks4^3*ks5*x01+6*ks3*ks4^4*ks5*x02)*_Z^5+(4*ks1^2*ks3^6*x02+8*ks1^2*ks3^5*ks4*x01-4*ks1^2*ks3^4*ks4^2*x02-4*ks1*ks2*ks3^5*x01^2+4*ks1*ks2*ks3^5*x02^2+16*ks1*ks2*ks3^4*ks4*x01*x02-8*ks1*ks2*ks3^3*ks4^2*x02^2-3*ks1*ks3^4*ks5*x01^2*x02+4*ks1*ks3^4*ks5*x02^3-2*ks1*ks3^3*ks4*ks5*x01^3+14*ks1*ks3^3*ks4*ks5*x01*x02^2+3*ks1*ks3^2*ks4^2*ks5*x01^2*x02-7*ks1*ks3^2*ks4^2*ks5*x02^3-ks1*ks4^4*ks5*x02^3-4*ks2^2*ks3^4*x01^2*x02-4*ks2^2*ks3^3*ks4*x01^3+4*ks2^2*ks3^3*ks4*x01*x02^2+4*ks2^2*ks3^2*ks4^2*x01^2*x02-3*ks2*ks3^3*ks5*x01^4-9*ks2*ks3^3*ks5*x01^2*x02^2-8*ks2*ks3^2*ks4*ks5*x01^3*x02+6*ks2*ks3^2*ks4*ks5*x01*x02^3+13*ks2*ks3*ks4^2*ks5*x01^2*x02^2+3*ks2*ks3*ks4^2*ks5*x02^4-2*ks2*ks4^3*ks5*x01*x02^3+12*ks3^9+24*ks3^7*ks4^2+12*ks3^5*ks4^4-7*ks3^2*ks5^2*x01^4*x02-7*ks3^2*ks5^2*x01^2*x02^3+6*ks3*ks4*ks5^2*x01^3*x02^2+6*ks3*ks4*ks5^2*x01*x02^4+ks4^2*ks5^2*x01^2*x02^3+ks4^2*ks5^2*x02^5+4*ks1*ks3^7+4*ks1*ks3^5*ks4^2-12*ks2*ks3^5*ks4*x01+12*ks2*ks3^4*ks4^2*x02+4*ks2*ks3^3*ks4^3*x01-4*ks2*ks3^2*ks4^4*x02+11*ks3^5*ks5*x01^2-42*ks3^4*ks4*ks5*x01*x02-17*ks3^3*ks4^2*ks5*x01^2+31*ks3^3*ks4^2*ks5*x02^2+38*ks3^2*ks4^3*ks5*x01*x02-21*ks3*ks4^4*ks5*x02^2)*_Z^4+(ks1^3*ks3^4*x01^2+ks1^3*ks3^4*x02^2+2*ks1^2*ks2*ks3^3*x01^2*x02+2*ks1^2*ks2*ks3^3*x02^3+2*ks1^2*ks3^2*ks5*x01^2*x02^2+2*ks1^2*ks3^2*ks5*x02^4+ks1*ks2^2*ks3^2*x01^2*x02^2+ks1*ks2^2*ks3^2*x02^4+2*ks1*ks2*ks3*ks5*x01^2*x02^3+2*ks1*ks2*ks3*ks5*x02^5+12*ks1*ks3^7*x02+28*ks1*ks3^6*ks4*x01-20*ks1*ks3^5*ks4^2*x02-4*ks1*ks3^4*ks4^3*x01+ks1*ks5^2*x01^2*x02^4+ks1*ks5^2*x02^6-5*ks2*ks3^6*x01^2+12*ks2*ks3^6*x02^2+50*ks2*ks3^5*ks4*x01*x02+11*ks2*ks3^4*ks4^2*x01^2-33*ks2*ks3^4*ks4^2*x02^2-14*ks2*ks3^3*ks4^3*x01*x02+3*ks2*ks3^2*ks4^4*x02^2-15*ks3^5*ks5*x01^2*x02+12*ks3^5*ks5*x02^3+2*ks3^4*ks4*ks5*x01^3+86*ks3^4*ks4*ks5*x01*x02^2+27*ks3^3*ks4^2*ks5*x01^2*x02-59*ks3^3*ks4^2*ks5*x02^3-44*ks3^2*ks4^3*ks5*x01*x02^2+15*ks3*ks4^4*ks5*x02^3+2*ks1^2*ks3^5*x02+2*ks1^2*ks3^4*ks4*x01-3*ks1*ks2*ks3^4*x01^2+2*ks1*ks2*ks3^4*x02^2+4*ks1*ks2*ks3^3*ks4*x01*x02+ks1*ks2*ks3^2*ks4^2*x02^2-2*ks1*ks3^3*ks5*x01^2*x02+2*ks1*ks3^3*ks5*x02^3+4*ks1*ks3*ks4^2*ks5*x02^3-4*ks2^2*ks3^3*x01^2*x02+4*ks2^2*ks3^2*ks4*x01*x02^2-9*ks2*ks3^2*ks5*x01^2*x02^2+8*ks2*ks3*ks4*ks5*x01*x02^3+ks2*ks4^2*ks5*x02^4+16*ks3^8+16*ks3^6*ks4^2-6*ks3*ks5^2*x01^2*x02^3+6*ks4*ks5^2*x01*x02^4+ks1*ks3^6+ks1*ks3^4*ks4^2-4*ks2*ks3^4*ks4*x01+4*ks2*ks3^3*ks4^2*x02+3*ks3^4*ks5*x01^2-12*ks3^3*ks4*ks5*x01*x02+9*ks3^2*ks4^2*ks5*x02^2)*_Z^3+(3*ks1^2*ks3^5*x01^2+3*ks1^2*ks3^5*x02^2+7*ks1*ks2*ks3^4*x01^2*x02+6*ks1*ks2*ks3^4*x02^3-ks1*ks2*ks3^2*ks4^2*x02^3+6*ks1*ks3^3*ks5*x01^2*x02^2+6*ks1*ks3^3*ks5*x02^4+4*ks1*ks3^2*ks4*ks5*x01*x02^3-4*ks1*ks3*ks4^2*ks5*x02^4+5*ks2^2*ks3^3*x01^2*x02^2+3*ks2^2*ks3^3*x02^4-2*ks2^2*ks3^2*ks4*x01*x02^3+11*ks2*ks3^2*ks5*x01^2*x02^3+6*ks2*ks3^2*ks5*x02^5-4*ks2*ks3*ks4*ks5*x01*x02^4-ks2*ks4^2*ks5*x02^5-4*ks3^8*x02+20*ks3^7*ks4*x01-36*ks3^6*ks4^2*x02-12*ks3^5*ks4^3*x01+7*ks3*ks5^2*x01^2*x02^4+3*ks3*ks5^2*x02^6-4*ks4*ks5^2*x01*x02^5+10*ks1*ks3^6*x02+10*ks1*ks3^5*ks4*x01-4*ks1*ks3^4*ks4^2*x02-3*ks2*ks3^5*x01^2+10*ks2*ks3^5*x02^2+22*ks2*ks3^4*ks4*x01*x02-13*ks2*ks3^3*ks4^2*x02^2-8*ks3^4*ks5*x01^2*x02+10*ks3^4*ks5*x02^3+38*ks3^3*ks4*ks5*x01*x02^2-24*ks3^2*ks4^2*ks5*x02^3+7*ks3^7+3*ks3^5*ks4^2)*_Z^2+(3*ks1*ks3^6*x01^2-4*ks1*ks3^6*x02^2-6*ks1*ks3^5*ks4*x01*x02+3*ks1*ks3^4*ks4^2*x02^2+5*ks2*ks3^5*x01^2*x02-4*ks2*ks3^5*x02^3-14*ks2*ks3^4*ks4*x01*x02^2+9*ks2*ks3^3*ks4^2*x02^3+7*ks3^4*ks5*x01^2*x02^2-4*ks3^4*ks5*x02^4-22*ks3^3*ks4*ks5*x01*x02^3+15*ks3^2*ks4^2*ks5*x02^4+ks1^2*ks3^4*x02^2+2*ks1*ks2*ks3^3*x02^3+2*ks1*ks3^2*ks5*x02^4+ks2^2*ks3^2*x02^4+2*ks2*ks3*ks5*x02^5-4*ks3^7*x02+8*ks3^6*ks4*x01-12*ks3^5*ks4^2*x02+ks5^2*x02^6+2*ks1*ks3^5*x02+2*ks2*ks3^4*x02^2+2*ks3^3*ks5*x02^3+ks3^6)*_Z-ks1^2*ks3^4*x02^3-2*ks1*ks2*ks3^3*x02^4-2*ks1*ks3^2*ks5*x02^5-ks2^2*ks3^2*x02^5-2*ks2*ks3*ks5*x02^6+ks3^7*x01^2-6*ks3^6*ks4*x01*x02+9*ks3^5*ks4^2*x02^2-ks5^2*x02^7-2*ks1*ks3^5*x02^2-2*ks2*ks3^4*x02^3-2*ks3^3*ks5*x02^4-ks3^6*x02)

where _Z stands in for the values that make the polynomial 0 -- the roots of the polynomial. You can break this up into parts that are the individual powers of _Z, get the 10 coefficient groups, and build a static roots() out of it.

Answer 1

[Bruno Luong]

It looks like you solve for pixel position of some distorded images, and you has a model for a distorsion.

If you approximate the inverse of the model, either symbolically or approximate by some sort of truncated taylor expansion of the inverse, might be you don't need to use FSOLVE at all.

[Yuk Ching Leung]

Yes, I am exactly working with distorted image and inversely calculating the undistorted position. You have a point where this problem can be solved differently using the inverse of the model. Thank you!