@maxburakho Your system of equations is linear and nonsingular. For such a system, Newton's method is guaranteed to converge to the exact solution in a single iteration. Newton's method involves solving a system of linear equations at each iteration. If the original equations are linear, that's the linear system that's solved anyway.

@Earl 
 

Download Operators.mw

Does that answer your question?

@lucaud The last thing in my already-posted Answer is the triple integral in cartesian coordinates for the volume. Is this not what you need? In set notation (which doesn't have much practical value in Maple), the volume can be described as

{(x,y,z)  in R^3 | -sqrt(1-y^2) <= z <= sqrt(1-y^2), -sqrt(1-x^2) <= y <= sqrt(1-x^2), -1 <= x <= 1}

OMG, I was wrong. Thanks for setting me straight. For positive b, or abs(b), the parabola's vertex and maximum is in the first quadrant.

@Kitonum The points G and H in the original diagram are part of the boundary of the convex hull, but they are not wanted by the OP.

For real b, maximize((b-x)*x, x) is b^2/4, but adding the domain restriction x >= 0 changes that! The maximum of ((abs(b)-x)*x) on x= 0..infinity is easily seen to be 0 because the parabola's vertex is in the second quadrant and the parabola passes through the origin. Unfortunately, Maple's maximize doesn't understand this:

maximize((abs(b)-x)*x, x= 0..infinity) assuming b::real;

returns unevaluated.

I don't fully understand your actual practical problem. In particular, I don't know if it has symbolic constants (such as the b above). If the constants have numeric values, I think that you can solve the problem with Optimization:-QPSolve.

Awesome! It hadn't occurred to me that morphing polygons would be so easy.

@vv I didn't mean to imply that extending to rational functions would involve simply changing A to rational functions. You need to multiply each row by the LCM (or just the product) of its denominators, make an estimate for the upper bound of the degree of the determinant of the new matrix (you can use the sum of the max of the degrees in each row), compute the determinant by the present process (with the number of imaging points determined by that degree estimate), and divide the result by the product of the LCMs.

@Rouben Rostamian  Yet the following shows how to using imaging and interpolation to compute the determinant of a 170 x 170 dense matrix of linear univariate polynomials with integer coefficients in 40 seconds:

restart:
n:= 170:
A:= LinearAlgebra:-RandomMatrix((n$2), generator= (()-> randpoly(x, degree= 1))):
Af:= (x::integer)-> rtable(eval(A, :-x= x), datatype= integer):
DET1:= CodeTools:-Usage(
   CurveFitting:-PolynomialInterpolation(
      [seq([k,LinearAlgebra:-Determinant(Af(k))], k= -iquo(n,2)..iquo(n,2)+1)],
      x
   )
):

memory used=0.82GiB, alloc change=48.03MiB, cpu time=39.88s, real time=40.29s, gc time=218.75ms

This can be easily distributed over multiple cores, like this (my present machine has 4 cores):

save DET1, `DET1.m`: #So we can compare results.
restart:
n:= 170:
A:= LinearAlgebra:-RandomMatrix((n$2), generator= (()-> randpoly(x, degree= 1))):
DET2:= CodeTools:-Usage(
   CurveFitting:-PolynomialInterpolation(
      [Grid:-Seq['tasksize'= 1](
         [k, LinearAlgebra:-Determinant(rtable(eval(A, x= k), datatype= integer))], 
         k= -iquo(n,2)..iquo(n,2)+1)
      ], x
   )
):

memory used=23.61MiB, alloc change=0 bytes, cpu time=1.02s, real time=18.02s, gc time=625.00ms

(That "memory used" and "cpu time" must only be measured for the master process, but the "real time" is real.)

read `DET1.m`:
evalb(DET1 = DET2);

     true

Extending this idea to A's entries being univariate rational functions with rational coefficients is not significantly more complicated.

This is essentially Rouben's method, but using only one line of code after the solve:

((F,B)-> eval(B, lhs~(F) =~ 0))(selectremove(evalb, sol_full));

Please post that code in plaintext so that I can copy-and-paste it. If you copy-and-paste directly from Maple to MaplePrimes, and you don't select "as image", that should do it.

@grzybs It's unfortunate that the conversion of a definite iterated integral to a different coordinate system requires several meticulous steps.

@Ronan Because of the different ways that mixtures of Pi and floats are treated in different versions of Maple, the Pi should be inside the brackets of the evalf. And you might as well put the 180 in there also. So,

evalf(FindAngle(p1,p2)*180/Pi);

@grzybs It looks like you forgot to copy-and-paste your new values, and simply retained Kitonum's values in your Reply above. Anyway, have a look at my worksheet below. The only thing that you need to change is the second-to-last command, which gives values to the constants.

So, did you find the above model useful or not? Some feedback would be appreciated.