applying vpa function for fsolve

%I need high precision while using fsolve function. My equations are badly scaled.

%I need to have like 100 digits after the decimal point. How can I perform this in the below codes?

F=@(x0,y0,z0)[(-4.71777*10^21)+(2.10263*(10^15)*x0)-(3.49196*(10^8)*(x0^2))+(28*x0^3)+(1.56321*(10^14)*y0)-(3.76035*(10^7)*x0*y0)-(1.16399*(10^8)*y0^2)+(28*x0*y0^2)+(1.11156*(10^15)*z0)-(2.6739*(10^8)*x0*z0)-(1.16399*(10^8)*z0^2)+(28*x0*z0^2)
    (-7.62057*10^20)+(1.56321*(10^14)*x0)-(1.88017*(10^7)*x0^2)+(1.16012*(10^15)*y0)-(2.32797*(10^8)*x0*y0)+(28*(x0^2)*y0)-(5.64052*(10^7)*y0^2)+(28*y0^3)+(1.7955*(10^14)*z0)-(2.6739*(10^8)*y0*z0)-(1.88017*(10^7)*z0^2)+(28*y0*z0^2)
    (-5.41882*(10^21))+(1.11156*(10^15)*x0)-(1.33695*(10^8)*x0^2)+(1.7955*(10^14)*y0)-(1.33695*(10^8)*y0^2)+(2.41161*(10^15)*z0)-(2.32797*(10^8)*x0*z0)+(28*(x0^2)*z0)-(3.76035*(10^7)*y0*z0)+(28*(y0^2)*z0)-(4.01085*(10^8)*z0^2)+(28*z0^3)
  ]

fun=@(u) F(u(1),u(2),u(3));

   x0=[4157000,671400,4774000]

[u1,fval1]=fsolve(fun,x0)

%actually fsolve function cannot find the solution therefore I want to try high precison, maybe it affects the solution.

Answer 1

[Walter Roberson]

You only have 6 digits of precision on your constants. You cannot meaningfully get out 100 digits of precision for the solution.

How close to 0 does F(u) need to be? For example, 25 digits gets F(u) to [-0.12e-1, -0.407e-4, -0.270e-2]

x0, y0, z0 have values that are primarily based on

roots([8730629105141324856992021631, -93154703113022319253101652905801, +5290130260748861049670439977885516000, -31683393121913492781344590435395215781250, +76812469157230892953736692218437414062500000, -92730743959203172134758219842091308593750000000, +57514411629616551846078000573281250000000000000000, 15761383693267952488406258764282226562500000000000000])

but of course carried out to more digits than the default.

[sermet]

dear Walter. I cannot understand what do you meant. Could you just show me the way that how can I apply vpa any digits I want before any calculation?

[Walter Roberson]

No. vpa is only used for symbolic calculations, but fsolve() is only or numeric calculations. Any extra digits you appear to get using vpa() will be lost when the symbolic number is converted to double precision that fsolve() works in.

If you want more digits you need to use the symbolic toolbox's numeric::fsolve

syms x0 y0 z0 real
F = sym( '[(-4.71777*10^21)+(2.10263*(10^15)*x0)-(3.49196*(10^8)*(x0^2))+(28*x0^3)+(1.56321*(10^14)*y0)-(3.76035*(10^7)*x0*y0)-(1.16399*(10^8)*y0^2)+(28*x0*y0^2)+(1.11156*(10^15)*z0)-(2.6739*(10^8)*x0*z0)-(1.16399*(10^8)*z0^2)+(28*x0*z0^2), (-7.62057*10^20)+(1.56321*(10^14)*x0)-(1.88017*(10^7)*x0^2)+(1.16012*(10^15)*y0)-(2.32797*(10^8)*x0*y0)+(28*(x0^2)*y0)-(5.64052*(10^7)*y0^2)+(28*y0^3)+(1.7955*(10^14)*z0)-(2.6739*(10^8)*y0*z0)-(1.88017*(10^7)*z0^2)+(28*y0*z0^2), (-5.41882*(10^21))+(1.11156*(10^15)*x0)-(1.33695*(10^8)*x0^2)+(1.7955*(10^14)*y0)-(1.33695*(10^8)*y0^2)+(2.41161*(10^15)*z0)-(2.32797*(10^8)*x0*z0)+(28*(x0^2)*z0)-(3.76035*(10^7)*y0*z0)+(28*(y0^2)*z0)-(4.01085*(10^8)*z0^2)+(28*z0^3)]' );

digits(50);
soln = feval(symengine, 'numeric::solve', F, [x0 = 4157000, y0 = 671400, z0 = 4774000])

Answer 2

[John D'Errico]

As Walter points out, you cannot get out more information than you put in. So if your constants are known only to 6 significant digits, asking for 100 digits is simply silly.

For example, when you use the number -4.71777*10^21 (a better way to write it is -4.71777e21), MATLAB actually uses the approximation

-4717769999999999797068.9946669153869152069091796875

Remember that MATLAB uses double precision for its numbers.

However, there is another problem here. FSOLVE does NOT work in 100 digits of precision. It is written to use the standard floating point, so double precision. Just wanting more accuracy will not suffice, nor will trying to use VPA on a tool that is not designed to work with it.

I imagine that you lack the symbolic TB, else you might have asked the question differently. However, lets see what happens when we perturb only one of the parameters in this problem.

syms x0 y0 z0

delta = sym('1.e-6');

E1 = -4.71777e21 + 2.10263e15*x0 - 3.49196e8*x0^2 + 28*x0^3 + 1.56321e14*y0 - ...
  3.76035e7*x0*y0 - 1.16399e8*y0^2 + 28*x0*y0^2 + 1.11156e15*z0 - 2.6739e8*x0*z0 - ...
  1.16399e8*z0^2 + 28*x0*z0^2;

E2 = -7.62057e20 + 1.56321e14*x0 - 1.88017e7*x0^2 + 1.16012e15*y0 - 2.32797e8*x0*y0 + ...
  28*x0^2*y0 - 5.64052e7*y0^2 + 28*y0^3 + 1.7955e14*z0 - 2.6739e8*y0*z0 - ...
  1.88017e7*z0^2 + 28*y0*z0^2;

E3 = -5.41882e21 + 1.11156e15*x0 - 1.33695e8*x0^2 + 1.7955e14*y0 - 1.33695e8*y0^2 + ...
  2.41161e15*z0 - 2.32797e8*x0*z0 + 28*x0^2*z0 - 3.76035e7*y0*z0 + 28*y0^2*z0 - ...
  4.01085e8*z0^2 + 28*z0^3;

res0 = solve(E1 == 0,E2 == 0,E3 == 0,x0,y0,z0);

% Perturbed problem
E1_p = -4.71777e21*(1 + sym('1.e-6')) + 2.10263e15*x0 - 3.49196e8*x0^2 + 28*x0^3 + 1.56321e14*y0 - ...
  3.76035e7*x0*y0 - 1.16399e8*y0^2 + 28*x0*y0^2 + 1.11156e15*z0 - 2.6739e8*x0*z0 - ...
  1.16399e8*z0^2 + 28*x0*z0^2;

res_p = solve(E1_p == 0,E2 == 0,E3 == 0,x0,y0,z0);

[res0.x0(1),res_p.x0(1)
res0.y0(1), res_p.y0(1)
res0.z0(1), res_p.z0(1)]
ans =
[ 4260085.7096479238399035675919043, 4265674.5969096452984634693062479]
[ 672704.55542203311984705958595163, 672651.74290433976042869910702819]
[  4849391.855485320782592481347078, 4846233.3214132306741776312904764]

So a perturbation of one part in 1 million to one of the constant terms has resulted in changes of a bit less than 1% in one of the solutions. The point is, you have specified these coefficients to only that relative accuracy. And worse, I'll bet that you don't really know those coefficients as well as you think. I'll bet that those coefficients are not truly known to the precision you show.

So sadly, asking for 100 digits of precision in the solution is silly in the extreme.

[Walter Roberson]

My tests suggest that if those constants are all exact, and if you bring them into the symbolic toolbox as being exact in base 10 (possible if you sym() a quoted string), then using 25 digits for calculation gives you an error in "x0" of about 0.012 which is too large for me to feel comfortable with. 30 digits might be more than you need though.