4
$\begingroup$

I have the following Mathematica code to solve an equation coming from isoparameteric finite element for quadratic tetrahedron element. I need to get the iteroplation fast and accurate. I only find articles by iterative method and I want to know whether it is possible to get an analytical solution:

sol = Solve[ ( 2(1 - xi - eta - zeta)^2 - (1 - xi - eta - zeta)) x1 + (2xi^2 - xi) x2 
             + (2eta^2 - eta) x3 + (2zeta^2 - zeta) x4 + 4(1 - xi - eta - zeta) xi x5 
             + 4xi eta x6 + 4(1 - xi - eta - zeta) eta x7 + 4eta zeta x8 
             + 4(1 - xi - eta - zeta) zeta x9 + 4xi zeta x10 == x  && 
            ( 2(1 - xi - eta - zeta)^2 - (1 - xi - eta - zeta)) y1 + (2xi^2 - xi) y2 
            + (2eta^2 - eta) y3 + (2zeta^2 - zeta) y4 + 4(1 - xi - eta - zeta) xi y5 
            + 4xi eta y6 + 4(1 - xi - eta - zeta) eta y7 + 4eta zeta y8 
            + 4(1 - xi - eta - zeta) zeta y9 + 4xi zeta y10 == y  &&  
            ( 2(1 - xi - eta - zeta)^2 - (1 - xi - eta - zeta))z1 + (2xi^2 - xi) z2
            + (2eta^2 - eta) z3 + (2zeta^2 - zeta) z4 + 4(1 - xi - eta - zeta) xi z5
            + 4xi eta z6 + 4(1 - xi - eta - zeta) eta z7 + 4eta zeta z8 
            + 4(1 - xi - eta - zeta) zeta z9 + 4xi zeta z10 == z, { xi, eta, zeta}];
xx =   xi /. sol;
yy =  eta /. sol;
zz = zeta /. sol;

by collecting the coefficient of xi, eta and zeta, the original equation can be transformed to the following:

sol = Solve[{  a1 xi^2 + a2 xi + a3 eta^2 + a4 eta + a5 zeta^2 + a6 zeta 
             + a7 xi eta + a8 eta zeta + a9 zeta xi == a0,
               b1 xi^2 + b2 xi + b3 eta^2 + b4 eta + b5 zeta^2 + b6 zeta 
             + b7 xi eta + b8 eta zeta + b9 zeta xi == b10,
               c1 xi^2 + c2 xi + c3 eta^2 + c4 eta + c5 zeta^2 + c6 zeta 
             + c7 xi eta + c8 eta zeta + c9 zeta xi == c10}, { xi, eta, zeta}];
xx =   xi /. sol;
yy =  eta /. sol;
zz = zeta /. sol;

where

{a1 b1 c1} = 2( {x1 y1 z1} + {x2 y2 z2} - 2{x5 y5 z5})

and so on. However, it seems to hang there for tens of minutes and it maybe impossible to get the analytical solution. I noticed that in real engineering, not all facets are curve facets. By this fact, we can suppose that 

2x5 = x1 + x2
2x7 = x1 + x3
2x9 = x4 + x1

and so on. Therefore, the above equation is simplified as follows:

    sol = Solve[{ a2 xi + a4 eta + a6 zeta + a7 xi eta + a8 eta zeta + a9 zeta xi == a10,
                  b2 xi + b4 eta + b6 zeta + b7 xi eta + b8 eta zeta + b9 zeta xi == b10, 
                  c2 xi + c4 eta + c6 zeta + c7 xi eta + c8 eta zeta + c9 zeta xi == c10}, 
                { xi, eta, zeta}];
xx =   xi /. sol;
yy =  eta /. sol;
zz = zeta /. sol;

Is there any good method to solve this equation?

|improve this question
$\endgroup$
  • $\begingroup$ Hi Verbeia, Thanks for help me to edit this item. Tang Laoya $\endgroup$ – Tang Laoya Jan 15 '13 at 5:32
  • 1
    $\begingroup$ You're welcome. To format code, just select the block of code and press the button that looks like { }. Or put four or more spaces at the beginning of each line. More info here. $\endgroup$ – Verbeia Jan 15 '13 at 5:37
18
$\begingroup$

One shouldn't expect too much working purely on the symbolic level with so many symbolic variables. More appropriate approach is to assume special values for several symbolic constants ai, bi, ci and to make a preliminary survey. Another (even more) important step is a distinction between generic and complete description of the solution sets. If the former is fully satisfactory you can use Solve, otherwise you should add to Solve the option MaxExtraConditions -> All or you should even better work with Reduce. For more complete discussion read What is the difference between Reduce and Solve?

As a hint what one can expect let's use ContourPlot3D and substitute special values for all symbolic constants, e.g. : 

ContourPlot3D[{ x + 2 y + 3 z + x y - y z + 2 x z == 1,
               -x + y + 3 z - x y + y z - 3 x z == -1, 
                x - 3 y + 2 z - 2 x y - y z + z x == 0 }, 
              {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, 
              ContourStyle -> {Green, Cyan, Orange}, MeshFunctions -> {#3 &}]

enter image description here

So we have three second order surfaces in three dimensional space. The solutions can be found where all three submanifolds intersect. One expects from 0 up to8 solutions. In this case for a concise output let's use Reduce :

Reduce[{ x + 2 y + 3 z + x y - y z + 2 x z == 1,
        -x + y + 3 z - x y + y z - 3 x z == -1, 
         x - 3 y + 2 z - 2 x y - y z + z x == 0 }, {x, y, z}]

 ( x == Root[-210 + 54 #1 + 182 #1^2 - 2 #1^3 - 22 #1^4 + 3 #1^5 &, 1] ||
   x == Root[-210 + 54 #1 + 182 #1^2 - 2 #1^3 - 22 #1^4 + 3 #1^5 &, 2] ||
   x == Root[-210 + 54 #1 + 182 #1^2 - 2 #1^3 - 22 #1^4 + 3 #1^5 &, 3] || 
   x == Root[-210 + 54 #1 + 182 #1^2 - 2 #1^3 - 22 #1^4 + 3 #1^5 &, 4] ||
   x == Root[-210 + 54 #1 + 182 #1^2 - 2 #1^3 - 22 #1^4 + 3 #1^5 &, 5]) &&
   y == 8 - 5 x - 4 x^2 + x^3 && z == 1/7 (-63 + 16 x + 47 x^2 + 4 x^3 - 3 x^4)

These are exact solutions represented in terms of the Root objects. If you prefer numeric values of the solutions you can do this : 

{x, y, z} /. Solve[{ x + 2 y + 3 z + x y - y z + 2 x z == 1,
                    -x + y + 3 z - x y + y z - 3 x z == -1, 
                     x - 3 y + 2 z - 2 x y - y z + z x == 0}, {x, y, z}] // N

 {{-1.87792, -3.33946, 1.2717}, {-1.5772, 2.01237, -0.796747},
  {0.985202, 0.147754, -0.0883909}, {4.70874, 0.170559, -0.396261},
  {5.09452, 10.9347, -36.2284}}

Now we can go further putting one symbolic constant in Reduce, e.g. c :

Reduce[{ x + 2 y + 3 z + x y - y z + 2 x z == 1,
        -x + y + 3 z - x y + y z - 3 x z == -1, 
         x - 3 y + 2 z - 2 x y - y z + z x == c }, {x, y, z}]

(c == 1 && x == 1 && y == 0 && z == 0) ||
(c == 1 && (x == 1/3 (5 - Sqrt[106]) || x == 1/3 (5 + Sqrt[106])) &&
           (y == 1/9 (9 + 8 x - Sqrt[486 + 45 x + 64 x^2]) || 
            y == 1/9 (9 + 8 x + Sqrt[486 + 45 x + 64 x^2])) && 
            z == 1/7 (-8 y - 3 x y))    ||
          ((x == Root[-630 - 117 c + 18 c^2 + (162 + 30 c - 3 c^2) #1 + (546 + 56 c) #1^2
                    + (-6 + 4 c) #1^3 + (-66 - 3 c) #1^4 + 9 #1^5 &, 1] || 
            x == Root[-630 - 117 c + 18 c^2 + (162 + 30 c - 3 c^2) #1 + (546 + 56 c) #1^2
                     + (-6 + 4 c) #1^3 + (-66 - 3 c) #1^4 + 9 #1^5 &, 2] || 
            x == Root[-630 - 117 c + 18 c^2 + (162 + 30 c - 3 c^2) #1 + (546 + 56 c) #1^2
                     + (-6 + 4 c) #1^3 + (-66 - 3 c) #1^4 + 9 #1^5 &, 3] || 
            x == Root[-630 - 117 c + 18 c^2 + (162 + 30 c - 3 c^2) #1 + (546 + 56 c) #1^2
                     + (-6 + 4 c) #1^3 + (-66 - 3 c) #1^4 + 9 #1^5 &, 4] || 
            x == Root[-630 - 117 c + 18 c^2 + (162 + 30 c - 3 c^2) #1 + (546 + 56 c) #1^2
                     + (-6 + 4 c) #1^3 + (-66 - 3 c) #1^4 + 9 #1^5 &, 5]) &&
          -1 + c != 0 && y == ((-1 + x) (24 + 3 c + 9 x + c x - 3 x^2))/(3 (-1 + c)) &&
                         z == 1/7 (1 - c - 8 y - 3 x y))

One can use Solve, but the output represented in terms of replacement rules would be substantially more involved. Sometimes it may apear that replacement rules cannot describe the solution set at all and then we would have to rely only on Reduce. Even for one symbolic constant the solution set appears to be quite involved and it is clear that one can make further progress only on a case by case basis. At this point one could use Manipulate and ContourPlot3D for all symbolic constants to explore the system with respect to more specific issues in order to reveal another interesting aspects. Here is the way to go taking a look at the link above for further references.

Edit

If we use only one symbolic parameter (variable) it is relatively easy to get exact symbolic 1-parameter solutions which are represented by segments of curves in 3-dimensional space {x,y,z}.
For a more compact notation let's introduce the following function :

sol[ a2_, a4_, a6_, a7_, a8_, a9_, a10_, b2_, b4_, b6_, b7_, b8_, b9_, b10_,
     c2_, c4_, c6_, c7_, c8_, c9_, c10_] := {x, y, z} /.
Solve[{ a2 x + a4 y + a6 z + a7 x y + a8 y z + a9 x z == a10, 
        b2 x + b4 y + b6 z + b7 x y + b8 y z + b9 z x == b10, 
        c2 x + c4 y + c6 z + c7 x y + c8 y z + c9 z x == c10}, {x, y, z}]

and now we are ready to explore possible solutions more extensively. Let's introduce functions of one variable, e.g. 

sl1[c_] := sol[1, 2, 3, 1, -1, 2, 1, -1, 1, 3, -1, 1, -3, -1, 1, -3, 2, -2, -1, 1, c]
sl2[d_] := sol[1, 2, 3, 1, -1, 2, 1, -1, 1, 3, -1, 1, -3d, -1, 1, -3, 2, -2, -1, 1, 0]
sl3[a_] := sol[1, 2, 3, 1, -1, 2, 1, -1, 1, 3, -1, 1, -3, -1, 1, -3, 2, -2a, -1, 1, 0]
sl4[b_] := sol[1, 2b, 3, 1, -1, 2, 1, -1, 1, 3, -1, 1, -3, -1, 1, -3, 2, -2, -1, 1, 0]

Having defined those functions we can demonstrate the solutions graphically by using ParametricPlot3D to visualize the dependance on different variables of sol.

GraphicsGrid[
    Map[ ParametricPlot3D[ #[c], {c, -50, 50}, Evaluated -> True, BoxRatios -> {1, 1, 1}, 
                           AxesLabel -> Table[Style[k, Bold, 24], {k, {x, y, z}}]]&,
         {{sl1, sl2}, {sl3, sl4}}, {2}]]

enter image description here

There are a few solutions for every argument therefore we have a few curves represented by different colors. The full symbolic solution sol is 21-dimensional submanifold in 24-dimensional space {x, y, z, a2, a4, ..., b2, b4,...c2,..., c10}. There is no way to get a complete symbolic representation of such a solution.

|improve this answer
$\endgroup$
  • $\begingroup$ Dear Dr. Artes, Thank you very much for your detailed explaination. I noticed that it is impossible or meanless to get the solution of all 'symbolic' equations. However, I am trying to solve this equation by iteration: 1) give an initial guess of z=z0; 2) sol=Solve[{a2*x+a4*y+a6*z+a7*xy+a8*yz+a9*zx==a10,b2*x+b4*y+b6*z+b7*xy+b8*yz+b9*zx==b10},{x,y}]; xx=x/.sol yy=y/.sol 3) update z by c2*x+c4*y+c6*z+c7*xy+c8*yz+c9*z*x==c10 4) if convergence return, else goto 1) Now the problem is that how to estimate the convergence of this iterative method? Any better iterative method? Thanks. $\endgroup$ – Tang Laoya Jan 16 '13 at 3:32
  • $\begingroup$ @TangLaoya You are welcome. I cannot see where you can go by iteration. In your simplified system you had three equations, thus you expected a few points in 3D space. Those points are isolated generically. Threfore I doubt that assuming a priori z = z0 would be good idea or I misundersood your approach. $\endgroup$ – Artes Jan 16 '13 at 3:43
  • $\begingroup$ Hi Artes, in the iteration I only solve {x,y} which is easy to give an expression by mathematica (contains z). The initial guess z0 is used to calculate x and y and then use third equation to update z, and so on. I don't know what's the convergence of this iteration or is there any better iterative method. Thanks. $\endgroup$ – Tang Laoya Jan 16 '13 at 3:47
  • $\begingroup$ @TangLaoya The method you have in mind cannot be very reliable unless you are interested in especially restricted cases. If I can I'll update tomorow my answer to include more snapshots of the solution space. As I wrote above a good probe for this task would be Manipulate and ContourPlot3D for symbolic parameters. $\endgroup$ – Artes Jan 16 '13 at 4:00
  • 1
    $\begingroup$ @TangLaoya I believe you should edit your question to make it more informative (including more details you've mentioned in the comments) and then we will see if we can do more. I am almost sure we can. At this point you could take a look at this post describing e.g. a numeric approach : mathematica.stackexchange.com/questions/16574/… $\endgroup$ – Artes Jan 16 '13 at 4:45
6
$\begingroup$

You can substitute some variables using replacement rules like this:

 Collect[{a1*u^2 + a2*v^2 + a3*w^2 + a4*x^2 + a5*y^2 + a6*z^2 == a7, 
    b1*u^2 + b2*v^2 + b3*w^2 + b4*x^2 + b5*y^2 + b6*z^2 == b7, 
    c1*u^2 + c2*v^2 + c3*w^2 + c4*x^2 + c5*y^2 + c6*z^2 == c7} /. {u ->
      x + y, v -> y + z, w -> z + x} // Expand, {x, y, z}]

(* {(a1 + a3 + a4) x^2 + (a1 + a2 + a5) y^2 + 
   2 a2 y z + (a2 + a3 + a6) z^2 + x (2 a1 y + 2 a3 z) == a7,
   (b1 + b3 + b4) x^2 + (b1 + b2 + b5) y^2 + 
   2 b2 y z + (b2 + b3 + b6) z^2 + x (2 b1 y + 2 b3 z) ==  b7,
  (c1 + c3 + c4) x^2 + (c1 + c2 + c5) y^2 + 
   2 c2 y z + (c2 + c3 + c6) z^2 + x (2 c1 y + 2 c3 z) == c7} *)

And then you will see that this is effectively a sixth-order polynomial, for which there aren't in general analytic solutions.

|improve this answer
$\endgroup$
  • $\begingroup$ Hi Verbeia, thanks for your kindly reply. Why do you say it is a sixth-order polynomial? I think it is second-order polynomial. Anything I have misunderstood? Tang Laoya $\endgroup$ – Tang Laoya Jan 15 '13 at 5:39
  • 1
    $\begingroup$ Well, I am not a mathematician, but you have three variables each raised to up to the second power, so that is more complex/higher-order than a single variable raised to the second power, which is what I'd call a second-order polynomial. Try substituting in numbers for the a, b and c coefficients and see if you get a result then, using NSolve. $\endgroup$ – Verbeia Jan 15 '13 at 5:45
  • $\begingroup$ Hi Verbeia, thank you very much for your kindly reply. Can we find a transform to elimate the terms xy, yz and zx in above equations (and don't increase x^1, y^1 and z^1 terms)? I know that for a invertible matrix A we can find P and Q to let PA*Q=D, where D is a diagonal matrix. Thanks. $\endgroup$ – Tang Laoya Jan 15 '13 at 5:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.