# equation solving – Plot is fast, solve is slow – what optimizations can be made?

I’m trying to solve a 3D geometric problem in terms of 3 variables (rotational angles), and I think the end is finally in sight – quite literally, I can see on a 3D plot where my solution exists. When I tell Mathematica to solve for it, however, it computes all night with no apparent progress.

Disclaimers: I come from a programming background and understand optimization, but I’m VERY new to Mathematica and have no idea what kind of optimizations need to be made. I’m hoping you can help point me in the right direction.

I’ve defined two functions:

1. rotateVector – this takes an axis (must be a unit vector), an angle, and a vector to rotate, and rotates the vector about the axis by the angle. It’s a simple matrix multiplication, so I think it’s probably relatively optimized already. For my specific problem, I’m only working in the reals, and using the built-in Mathematica rotation seemed to add more complexity than I needed. If that was a dumb mistake, I’m happy to correct.
``````rotateVector(axis_, angle_, vector_) := Module(
{rotationMatrix, cosAngle, cosAngleC, sinAngle, ux, uy, uz},
ux = axis((1));
uy = axis((2));
uz = axis((3));
cosAngle = Cos(angle);
cosAngleC = 1 - cosAngle;
sinAngle = Sin(angle);
rotationMatrix = {{cosAngle + ux*ux*cosAngleC,
ux*uy*cosAngleC - uz*sinAngle,
ux*uz*cosAngleC + uy*sinAngle}, {uy*ux*cosAngleC + uz*sinAngle,
cosAngle + uy*uy*cosAngleC,
uy*uz*cosAngleC - ux*sinAngle}, {uz*ux*cosAngleC - uy*sinAngle,
uz*uy*cosAngleC + ux*sinAngle, cosAngle + uz*uz*cosAngleC}};
rotationMatrix.vector
)
``````
1. calcVectors2 (as in version 2) – this takes 3 arguments: an “outward angle” connected to the first transform, a “bend angle” connected to the 2nd, and a “crease angle”. It used to take a fourth, “crease amount,” but I found a way to solve for that.
``````calcVectors2(outAngle_, bendAngle_, creaseAngle_) := Module(
{p1, p1Out, p1Bend, crease, creaseOut, creaseBend, creaseMirror,
creaseMirrorProj, creaseOrth, p1Proj, p1Orth, creaseAmount,
afterCrease, afterCreaseMirror},
p1 = {0, 0, -1};
p1Out = rotateVector({1, 0, 0}, outAngle, p1);
p1Bend = rotateVector({0, 1, 0}, bendAngle, p1Out);
crease = rotateVector({0, 1, 0}, creaseAngle, p1);
creaseOut = rotateVector({1, 0, 0}, outAngle, crease);
creaseBend = rotateVector({0, 1, 0}, bendAngle, creaseOut);
creaseMirror = {-creaseBend((1)), creaseBend((2)),
creaseBend((3))};
creaseMirrorProj = Projection(creaseMirror, creaseBend);
creaseOrth = creaseMirror - creaseMirrorProj;
p1Proj = Projection(p1Bend, creaseBend);
p1Orth = p1Bend - p1Proj;
creaseAmount = -VectorAngle(p1Orth, creaseOrth);
afterCrease = rotateVector(creaseBend, creaseAmount, p1Bend);
afterCreaseMirror = {-afterCrease((1)), afterCrease((2)),
afterCrease((3))};
VectorAngle(afterCreaseMirror, creaseBend)
)
``````

If it helps to understand my problem: I take a unit vector pointing straight down in Z, called “p1.” I also generate a “crease vector” in the XZ plane (lower left quadrant), at angle “creaseAngle” to p1. I rotate both vectors about the X axis by “outward angle” (in the 0 to -Pi/2 range), then rotate about the Y axis by “bend angle” (0 to -Pi/4). Finally, I rotate the resultant “p1Bend” about the “creaseBend” by an amount I calculate to bring it into a plane defined between the “creaseBend” vector, and a mirror of that vector over the YZ plane. The final constraint: mirror the resultant vector that started as p1 over the YZ plane, and compare the angle between it and the “creaseBend” vector – I want it to be 0.

When I put in a fixed value for outAngle (-Pi/4), and “x” and “y” for my “bendAngle” and “creaseAngle”, I can see a plot where I can SEE the solution line – shaped roughly parabola-ish – this only takes a couple minutes or so to plot:

`Plot3D(calcVectors2(-Pi/4, x, y), {x, -((Pi)/3), (Pi)/ 12}, {y, -((Pi)/12), (7 (Pi))/12}, PlotPoints -> 70, MaxRecursion -> 7)` When I tell it to solve, however… it just chugs and chugs and chugs and nothing seems to happen; I gave up after 8 hours. I tried to simplify the problem, and assign a fixed value (-Pi/8) for x so I could see just a slice, and again the 0 is very visually apparent:

`Plot(calcVectors2(-Pi/4, -Pi/8, x), {x, -Pi/12, 7 Pi/12})` However, when I try to solve for `calcVectors2(-Pi/4, -Pi/8, x) == 0`, it just hangs again. Admittedly, I haven’t tried running this for 8 hours yet, but I figure I have some TERRIBLY unoptimized code, and I was wondering where I should look for speedups. I would love to be able to find a general solution for all 3 arguments to calcVectors2, but in terms of practical usage, the first 2 arguments are chosen for a given application. I guess what I’m ultimately after is a formula for the 3rd argument, using the first 2 arguments as input.

Where can I optimize so it doesn’t take forever to solve?

Posted on Categories Articles