Inaccurate answer for a system of equations

When I solve the following system of 6 equations by NSolve, it gives me an answer set; however this answer won't make the 2nd equation equal to zero (i.e. when I put the answers into

-π/6*dd^3 + π*(10^-3*L) ... == 0

it gives me a residual which is from the same order of the first term (-π/6*dd^3) Please let me know how I should change the code to get the exact solution of system.

Many thanks,

Clear[L, db, dd, dp, vdbi, vpbi, mp, md, sigma, mu, rod, rop, vpai, m, h, delt, H, vlam, c1, c2, c3, c4, θ]
ClearAll

dd = 2.83*10^-3;
dp = 2*10^-3;
vdbi = 0.82;
vpbi = 9;
mp = 4.4*10^-6;
md = 1.18*10^-5;
sigma = 0.0734;
rod = 1000;
rop = 1050;
mu = 0.005;
vpai = 6.71;
m = 3;

NSolve[{-π/12 (rod*dd^3*vdbi^2 + rop*dp^3*vpbi^2) + π*sigma*(10^-3*L) (2dp+ (10^-3*2 L)*Sqrt[1 - θ^2]) + π^2*(10^-4*2 h)*(dp + (10^-3*2*L)*Sqrt[1 - θ^2]) sigma + π*mu*(10^-3*L)*(2*dp + (10^-3*2*L)*Sqrt[1 - θ^2])/2*64/(h*10^-4)*(10^-3*L)*θ/vpai + π/12*rop*dp^3*vpai^2 + π/8*dp^2*rop*.5*(vpai*dd)^2/(10^-3*L*θ) == 0,
-π/6*dd^3 + π*(10^-3*L)/2*(2 dp + (10^-3*2 L)*Sqrt[1 - θ^2])*(10^-4*h) + π^2*(10^-4*h)^2 (dp + (10^-3*2 L)*Sqrt[1 - θ^2]) == 0,
vlam*(10^-3*L)*θ/vpai - 4 π*(10^-4*h) == 0,
6*(10^12*c1)/m*(10^-4*h) + 2 (c2*10^8) == 0 ,
(10^12*c1)*((10^-4*h/m)^3 - 3*10^-12*h^3/m) + (10^8*c2)*10^-8*h^2/m (1/m - 2) + vpai == 0,
(10^12*c1)*10^-12*h^3/4 + (10^8*c2)*10^-8*h^2/3 - ((10^12*1.5 c1)*10^-4*h + (10^8*c2))*10^-8*h^2 + vpai + vlam == 0,
h > 0,
L > 0,
θ > 0,
vlam > 0
}, {L, h, vlam, c1, c2, θ}, Reals]
• This is a numerical solution, you don't necessarily get exactly zero. You should decide that what is zero? ($10^{-3}$,$10^{-10}$ or $10^{-50}$), it depends on the problem. Jun 13 '15 at 22:50
• @Mahdi Thanks for your comment. But, for instance the 2nd equation consists of 3 terms and order of magnitude for zero that NSolve gives me is equal to the order of magnitude of each 3 term (for example, the first term, π/6*dd^3 , is equal to 1.18e-8 and what I got instead of zero from NSolve is equal to 2.27e-8 that indicates my answers can't be accurate at all. Do you know how I can enhance it? Thanks again Jun 13 '15 at 23:51
• I get { } when I run your code. Please check whether your code in the question is correct. Also, what is Clear All? Jun 14 '15 at 4:15
• @bbgodfrey Thanks for your comment. I run the code and it gives the following answer: {{vlam -> 2.8975, c2 -> -1.44137, c1 -> 0.426397, θ -> 0.995915, L -> 9.87748, h -> 3.38033}} However, I added the first line to clear all existing values. I appreciate if you can help me with this problem. Jun 14 '15 at 4:54
• I get the same as bbgodfrey. What version are you using? I'm using V10.1. Jun 14 '15 at 5:01

This seems a thorny and or buggy system. I dug out a V10.0.0 and reproduced the OP's result. I got somewhere with it on V10.1. Some changes caused the kernel to crash. It may also be that there are no real solutions to the system. The following code was carried out on V10.1.

First some definitions to make it easier to deal with exploring the issues:

sys = {-π/12 (rod*dd^3*vdbi^2 + rop*dp^3*vpbi^2) + π*
sigma*(10^-3*L) (2 dp + (10^-3*2 L)*
Sqrt[1 - θ^2]) + π^2*(10^-4*2 h)*(dp + (10^-3*2*L)*Sqrt[1 - θ^2]) sigma + π*
mu*(10^-3*L)*(2*dp + (10^-3*2*L)*Sqrt[1 - θ^2])/
2*64/(h*10^-4)*(10^-3*L)*θ/vpai + π/12*rop*dp^3*
vpai^2 + π/8*dp^2*rop*.5*(vpai*dd)^2/(10^-3*L*θ) == 0,
-π/6*dd^3 + π*(10^-3*L)/
2*(2 dp + (10^-3*2 L)*Sqrt[1 - θ^2])*(10^-4*
h) + π^2*(10^-4*h)^2 (dp + (10^-3*2 L)*Sqrt[1 - θ^2]) == 0,
vlam*(10^-3*L)*θ/vpai - 4 π*(10^-4*h) == 0,
6*(10^12*c1)/m*(10^-4*h) + 2 (c2*10^8) == 0,
(10^12*c1)*((10^-4*h/m)^3 - 3*10^-12*h^3/m) + (10^8*c2)*10^-8*
h^2/m (1/m - 2) + vpai == 0,
(10^12*c1)*10^-12*h^3/4 + (10^8*c2)*10^-8*
h^2/3 - ((10^12*1.5 c1)*10^-4*h + (10^8*c2))*10^-8*h^2 + vpai +
vlam == 0,
h > 0, L > 0, θ > 0, vlam > 0};
vars = {L, h, vlam, c1, c2, θ};

sysrat = Rationalize[sys, 0];          (* system with exact coefficients *)
eqnrat = Cases[sysrat, _Equal];        (* equations *)
resrat = eqnrat /. Equal -> Subtract;  (* residuals *)

Rationalizing the system allows us to play with the precision.

NSolve[sys, vars, Reals]
sol20 = NSolve[sysrat, vars, Reals, WorkingPrecision -> 20]
NSolve[sysrat, vars, Reals, WorkingPrecision -> 30]
(*
{}
{{L -> 9.877482896535161094, h -> 3.380334356733960619,
vlam -> 2.8975000000000000000, c1 -> 0.4263974715325724480,
c2 -> -1.4413660226460455712, \[Theta] -> 0.9959150897330469926}}
{}
*)

While at MachinePrecision, NSolve indicates there is no solution, at a precision of 20, it reproduces the OP's result from V10.0.0, including the large first residual:

resrat /. sol20
(*
{{0.000086692384737759982, 2.27342360012992162*10^-8, 0.*10^-22,
0.*10^-11, 0.*10^-18, 0.*10^-18}}
*)

On the other hand, at higher precisions, NSolve goes back to reporting no solution. Further, each of the following causes the kernel to crash:

NSolve[sysrat, vars, Reals]
NSolve[eqnrat, vars, Reals]
NSolve[eqnrat, vars]

Crashing the kernel seems like a bug to me, but I have not investigated further.

A common approach to try when equations do not behave tamely is to minimize the sum of the squares of the residuals. I had to give NMinimize some help. By itself (using the defaults), it was far from the best solution, producing a sum of squares equal to 983555.. (The results here are the same on V10.0.0.)

{resSS, solm} =
NMinimize[{resrat^2 // Total, DeleteCases[sysrat, _Equal]}, vars,
"InitialPoints" ->
Join[vars /. sol20,
vars + 10^-1 IdentityMatrix[Length[vars]] /. First@sol20]}]
(*
{7.51557*10^-9,
{L -> 9.87748, h -> 3.38033, vlam -> 2.8975,
c1 -> 0.426397, c2 -> -1.44137, θ -> 0.995915}}
*)

This is the same solution as NSolve. In other words, it cannot improve upon the NSolve result. So I'm led to guess that there are no real solutions to the OP's system.

Note on V10.0.0

If we relax the constraint, we get two more real solutions, which may actually be solutions, using one of the commands that crash the kernel in V10.1:

NSolve[eqnrat, vars, Reals]
(*
{{vlam -> 2.8975, c2 -> -1.44137, c1 -> 0.426397,
θ -> 0.995915, L -> 9.87748, h -> 3.38033},
{vlam -> 2.8975, c2 -> -0.659392, c1 -> 0.131938,
θ -> -0.980501, L -> -14.8332, h -> 4.99775},
{vlam -> 2.8975, c2 -> -0.0467831, c1 -> 0.00249337,
θ -> -0.999097, L -> -54.6517, h -> 18.763}}
*)

resrat /. % // Chop[#, 10^-16] &
(*
{{0.0000866924, 2.27342*10^-8, 0, -5.36442*10^-7, 0, 0},
{0, 0, 0, -2.38419*10^-7, 0, 0},
{0, 0, 0, -1.49011612`*^-8, 0, 0}}
*)
• Thanks for your helpful answer. In the answer set I got from NSolve, for 5 variables (L, h, vlam, c1, c2) it returns accurate values and only the value found for θ is not accurate. In this case, whether I can trust on those 5 values rather than θ that I got? (or because θ has a large error, they are also "definitely" unacceptable?) Thanks Jun 15 '15 at 12:32
• @Nick You're welcome. But I'm confused about your remark. How do you know that only θ is off? What should it be? (Do you already know the solution and you're just checking NSolve? Including the solution in the question would have been helpful in analyzing the problem, if it were known.) Jun 15 '15 at 15:14
• oh yes. These equations are related to a physical system and since we've performed some experiments on that as well, we can compare the results. Interestingly, the values we get for L, h, vlam, c1, c2 are close to what we anticipate based on experiments, while for θ it's not. (we anticipate θ to be around 0.93). Thanks again Jun 15 '15 at 18:18