calculation speed

• To: mathgroup at smc.vnet.net
• Subject: [mg71880] calculation speed
• From: "Jean-Paul" <Jean-Paul.VINCENT-3 at etudiants.ensam.fr>
• Date: Sat, 2 Dec 2006 05:10:46 -0500 (EST)

```Dear Group,
In order to solve a non linear system (6*6), i calculate with
Mathematica the Jacobian Matrix of the system.

The resolution's algorithm needs to calculate the inverse of this
Matrix.

Expressions i handle are very complex. So it takes 1.5-1.9 sec to
Mathematica too calculate it.

How can i reduce this time?

Here is my program:

Remove["Global`*"];
Off[Part::"pspec"]
Off[General::"spell"]
Needs["Optimization`UnconstrainedProblems`"]

(*Introduction des matrices de rotation et des erreurs*)
z11 = 13;
z22 = 21;
Î± = Pi/7.5;
(*Angle en les axes de rotation*)
Ï?0 = Pi/2;
coeftaille = 1;
{Rp, Rg, B,
A1, B1, C1, A2, B2, C2} = {50*coeftaille, 100*coeftaille, \
0.11*coeftaille, 0.02, 2.5, 0.05, 0.02, 2.5, 0.05};
omegamaxi = B/((2*Rg + Rp)/3*Sin[Î´b1]);
Î´1 = ArcTan[Sin[Ï?0]/(z22/z11 + Cos[Ï?0])];
Î´2 = Ï?0 - Î´1;
Î´b1 = ArcSin[Sin[Î´1]*Cos[Î±]];
Î´b2 = ArcSin[Sin[Î´2]*Cos[Î±]];
a1 = Sin[Î´b1];
b1 = Cos[Î´b1];
a2 = Sin[Î´b2];
b2 = Cos[Î´b2];
erreurpasp = {0, 0.0002, -0.0002, 0, -0.0001, +0.0001,
0.0005, -0.0005, 0, -0.000002, 0.000002, 0, 0};
erreurpasr = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0};

(*Erreur angulaire entre ces axes*)
Ï? = 0.00025;
(*Angle rÃ©el entre les axes*)
Ï? = Ï?0 + Ï?;
(*autres erreurs*)
Î?1 = 10^-4;
u35 = 0;
v35 = 0;
w35 = 0;
A11 = 0;
A22 = 0;
A33 = 0;
Î?2 = 0;
u46 = 0;
v46 = 0;
w46 = 0;
npign = -2*Pi/z11*(dent - 1) + Sum[erreurpasp[[i]], {i, 1, dent}];
nroue = 2*Pi/z22*(dent - 1) + Sum[erreurpasr[[i]], {i, 1, dent}];
lamma1 = ArcCos[Cos[Î´1]/Cos[Î´b1]];
Î´a1 = Î´1 + Î½a1;
Î½a1 = Ka1*num;
num = 2*Sin[Î´1]/z11;
Ka1 = Kh - Kfe1;
Kh = 2;
Kfe1 = Kfe1maxi - cKfe1;
Kfe1maxi = (Î´1 - Î´b1)/num;
cKfe1 = 0.05;
T10S1 = ArcCos[Cos[Î´a1]/Cos[Î´b1]];
ISS1 = T10S1 - lamma1;
I01S1 = ArcCos[(Cos[ISS1] - Cos[Î´1]*Cos[Î´a1])/(Sin[Î´1]*Sin[Î´a1])];
Î¸S1 = 1/a1*ArcCos[Cos[Î´a1]/b1];
p2 = ArcTan[(a1*Sin[Î¸S1]*Cos[a1*Î¸S1] -
Cos[Î¸S1]*Sin[a1*Î¸S1])/(a1*Cos[
Î¸S1]*Cos[a1*Î¸S1] + Sin[Î¸S1]*Sin[a1*Î¸S1])];
Î¦1mini = -I01S1 - p2;

lamma2 = ArcCos[Cos[Î´2]/Cos[Î´b2]];
Î´a2 = Î´2 + Î½a2;
Î½a2 = Ka2*num;
Ka2 = Kh - Kfe2;
Kfe2 = Kfe2maxi - cKfe2;
Kfe2maxi = (Î´2 - Î´b2)/num;
cKfe2 = 0.05;
T20S2 = ArcCos[Cos[Î´a2]/Cos[Î´b2]];
ISS2 = T20S2 - lamma2;
O1SS2 = ArcCos[Cos[ISS2]*Cos[Î´1] + Sin[ISS2]*Sin[Î´1]*Cos[Pi/2 - Î±]];
Î¸S2 = 1/a1*ArcCos[Cos[O1SS2]/b1];
p1 = ArcTan[(a1*Sin[Î¸S2]*Cos[a1*Î¸S2] -
Cos[Î¸S2]*Sin[a1*Î¸S2])/(a1*Cos[Î¸S2]*Cos[a1*Î¸S2] + Sin[
Î¸S2]*Sin[a1*Î¸S2])];
I01S2 = ArcCos[(Cos[Î´a2] -
Cos[Ï?0]*Cos[O1SS2])/(Sin[Ï?0]*Sin[O1SS2])];
Î¦1maxi = I01S2 - p1;

In[66]:=

(*Pignon*)
x1[u1_,v1_]=v1*(b1*Cos[a1*u1]);
y1[u1_,v1_]=v1*(a1*Cos[u1]*Cos[a1*u1]+Sin[u1]*Sin[a1*u1]);
z1[u1_,v1_]=v1*(a1*Sin[u1]*Cos[a1*u1]-Cos[u1]*Sin[a1*u1]);
s1[u1_,v1_]={x1[u1,v1],y1[u1,v1],z1[u1,v1]}//FullSimplify;

(*Roue*)
x2[u2_,v2_]=v2*(b2*Cos[a2*u2]);
y2[u2_,v2_]=v2*(a2*Cos[u2]*Cos[a2*u2]+Sin[u2]*Sin[a2*u2]);
z2[u2_,v2_]=v2*(a2*Sin[u2]*Cos[a2*u2]-Cos[u2]*Sin[a2*u2]);
s2[u2_,v2_]={x2[u2,v2],y2[u2,v2],z2[u2,v2]}//FullSimplify;

In[74]:=

(*Ajout du bombÃ© sur le pignon*)

Simplify;

In[76]:=
(*calcul de la normale pour cette gÃ©ometrie*)
sb1u1[u1_,v1_]=D[sb1[u1,v1],u1]//Simplify;
sb1v1[u1_,v1_]=D[sb1[u1,v1],v1]//Simplify;
nb1[u1_,v1_]={sb1u1[u1,v1][[2]]*sb1v1[u1,v1][[3]]-sb1u1[u1,v1][[3]]*sb1v1[u1,
v1][[2]],sb1u1[u1,v1][[3]]*sb1v1[u1,

v1][[1]]-sb1u1[u1,v1][[1]]*sb1v1[u1,v1][[3]],sb1u1[u1,v1][[1]]*sb1v1[u1,\
v1][[2]]-sb1u1[u1,v1][[2]]*sb1v1[u1,v1][[1]]};
normnb1[u1_,v1_]=Sqrt[nb1[u1,v1][[1]]*nb1[u1,v1][[1]]+nb1[u1,v1][[2]]*nb1[
u1,v1][[2]]+nb1[u1,v1][[3]]*nb1[u1,v1][[3]]];
nb1u[u1_,v1_]=nb1[u1,v1]/normnb1[u1,v1];

In[81]:=
(*Calcul de la fonction dÃ©formation*)

(*pignon*)
w1[u1_,v1_]=u1^2+u1*v1*v1^2//Simplify;

(*roue*)
w2[u2_,v2_]=u2^2+u2*v2+v2^2//Simplify;

In[83]:=
(*Calcul de la normale pour la roue*)
s2u2[u2_,v2_]=D[s2[u2,v2],u2]//Simplify;
s2v2[u2_,v2_]=D[s2[u2,v2],v2]//Simplify;
n2[u2_,v2_]={s2u2[u2,v2][[2]]*s2v2[u2,v2][[3]]-s2u2[u2,v2][[3]]*s2v2[u2,v2][[
2]],s2u2[u2,v2][[3]]*s2v2[u2,v2][[1]]-s2u2[u2,
v2][[1]]*s2v2[
u2,v2][[3]],s2u2[u2,v2][[1]]*s2v2[u2,v2][[2]]-
s2u2[u2,v2][[2]]*s2v2[u2,v2][[1]]};
normn2[u2_,v2_]=Sqrt[n2[u2,v2][[1]]*n2[u2,v2][[1]]+n2[
u2,v2][[2]]*n2[u2,v2][[2]]+n2[u2,v2][[3]]*n2[u2,v2][[3]]];
n2u[u2_,v2_]=n2[u2,v2]/normn2[u2,v2];

In[88]:=
(*GÃ©ometries avec dÃ©fauts de forme dans les repÃ¨res R7 et R8*)

(*Pignon*)

xp1[u1_,v1_]=sb1[u1,v1][[1]]+nb1u[u1,v1][[1]]*w1[u1,v1];
yp1[u1_,v1_]=sb1[u1,v1][[2]]+nb1u[u1,v1][[2]]*w1[u1,v1];
zp1[u1_,v1_]=sb1[u1,v1][[3]]+nb1u[u1,v1][[3]]*w1[u1,v1];
sp1[u1_,v1_]={xp1[u1,v1],yp1[u1,v1],zp1[u1,v1]};

(*Roue*)
xr2[u2_,v2_]=s2[u2,v2][[1]]+n2u[u2,v2][[1]]*w2[u2,v2];
yr2[u2_,v2_]=s2[u2,v2][[2]]+n2u[u2,v2][[2]]*w2[u2,v2];
zr2[u2_,v2_]=s2[u2,v2][[3]]+n2u[u2,v2][[3]]*w2[u2,v2];
sr2[u2_,v2_]={xr2[u2,v2],yr2[u2,v2],zr2[u2,v2]};

(*Pignon*)
(*Matrice de rotation pour passer de la premiÃ¨re Ã  la Dent - iÃ¨me
dent*)
m57 = {{1, 0, 0, 0}, {0, Cos[npign], -Sin[npign], 0}, {0, Sin[
npign], Cos[npign], 0}, {0, 0, 0, 1}};
(*Matrice de rotation autour de y pour simuler un problÃ¨me
d'alignement*)
m35 = {{Cos[Î?1], 0, -Sin[Î?1], u35}, {0, 1, 0, v35},
{Sin[Î?1], 0, Cos[Î?1], w35}, {0, 0, 0, 1}};
(*translation des sommets des cÃ´nes*)
m01 = {{1, 0, 0, A11}, {0, 1, 0, A22}, {0, 0, 1, A33}, {0, 0, 0, 1}};
(*rotation de Î¦1 autour de x*)
m13 = {{1, 0, 0, 0}, {
0, Cos[Î¦1], -Sin[Î¦1], 0}, {0, Sin[Î¦1], Cos[Î¦1], 0}, {0, 0, 0, 1}};
(*Matrice globale*)
mpign = m01.m13.m35.m57;

(*Roue*)
(*Matrice de rotation pour passer de la premiÃ¨re Ã  la Dent - iÃ¨me
dent*)
m68 = {{1, 0, 0, 0}, {0, Cos[nroue], -Sin[nroue], 0}, {0, Sin[nroue], \
Cos[nroue], 0}, {0, 0, 0, 1}};
(*Matrice de rotation autour de y pour simuler un problÃ¨me
d'alignement*)
m46 = {{Cos[Î?2], 0, -Sin[Î?2], u46}, {0, 1, 0, v46},
{Sin[Î?2], 0, Cos[Î?2], w46}, {0, 0, 0, 1}};
(*matrice de rotation pour mettre en place l'axe du pignon avec celui
du \
bati*)
m02 = {{Cos[Ï?], -Sin[Ï?], 0, 0}, {Sin[Ï?], Cos[Ï?], 0, 0}, {0, 0, 1,
0}, {0, 0, \
0, 1}};
(*rotation de Î¦2 autour de x*)
m24 = {{1, 0, 0, 0}, {0, Cos[Î¦2], -Sin[Î¦2], 0}, {0, Sin[Î¦2], Cos[
Î¦2], 0}, {0, 0, 0, 1}};
(*Matrice globale*)
mroue = m02.m24.m46.m68;

(*Pignon*)
w1u1[u1_, v1_] = D[w1[u1, v1], u1] // Simplify;
nb1uu1[u1_, v1_] = D[nb1u[u1, v1], u1];
sp1u1[u1_,
v1_] = sb1u1[u1, v1] + w1u1[u1, v1]*nb1u[
u1, v1] + w1[u1, v1]*nb1uu1[u1, v1];
w1v1[u1_, v1_] = D[w1[u1, v1], v1] // Simplify;
nb1uv1[u1_, v1_] = D[nb1u[u1, v1], v1];
sp1v1[u1_,
v1_] = sb1v1[u1, v1] + w1v1[u1,
v1]*nb1u[u1, v1] + w1[u1, v1]*nb1uv1[u1, v1];
np1d[u1_, v1_] = {sp1u1[u1, v1][[2]]*sp1v1[u1, v1][[3]] - sp1u1[u1,
v1][[3]]*sp1v1[u1, v1][[2]],
sp1u1[u1, v1][[3]]*sp1v1[u1, v1][[1]] - sp1u1[u1,
v1][[1]]*sp1v1[u1,
v1][[3]], sp1u1[u1, v1][[1]]*sp1v1[u1,
v1][[2]] - sp1u1[u1, v1][[2]]*sp1v1[u1, v1][[1]]};

(*roue*)
w2u2[u2_, v2_] = D[w2[u2, v2], u2] // Simplify;
n2uu2[u2_, v2_] = D[n2u[u2, v2], u2];
sr2u2[u2_, v2_] = s2u2[u2, v2] + w2u2[u2, v2]*n2u[u2, v2] + w2[u2,
v2]*n2uu2[u2, v2];
w2v2[u2_, v2_] = D[w2[u2, v2], v2] // Simplify;
n2uv2[u2_, v2_] = D[n2u[u2, v2], v2];
sr2v2[u2_, v2_] = s2v2[u2, v2] + w2v2[u2, v2]*n2u[u2, v2] + w2[u2, v2]*
n2uv2[u2, v2];
nr2d[u2_, v2_] = {sr2u2[u2,
v2][[2]]*sr2v2[u2, v2][[3]] -
sr2u2[u2, v2][[3]]*sr2v2[u2, v2][[2]], sr2u2[u2, v2][[
3]]*sr2v2[u2, v2][[1]] - sr2u2[
u2, v2][[1]]*sr2v2[u2, v2][[
3]], sr2u2[u2, v2][[1]]*sr2v2[
u2, v2][[2]] - sr2u2[u2, v2][[2]]*sr2v2[u2, v2][[1]]};

(* Ecriture de la gÃ©omÃ©trie et des normales dans le repÃ¨re global
pour \
(*Pignon*)
tablepign = mpign.Flatten[{sp1[u1, v1], 1}];
sp11[u1_, v1_, dent_, Î¦1_] = {tablepign[[1]], tablepign[[2]],
tablepign[[3]]};

tableNpign = mpign.Flatten[{np1d[u1, v1], 1}];
np11[u1_, v1_, dent_, Î¦1_] = {tableNpign[[1]], tableNpign[[2]],
tableNpign[[
3]]};

(*Roue*)
tableroue = mroue.Flatten[{sr2[u2, v2], 1}];
sr22[u2_, v2_, dent_, Î¦2_] = {tableroue[[1]], tableroue[[2]],
tableroue[[3]]};
tableNroue = mroue.Flatten[{nr2d[u2, v2], 1}];
nr22[u2_, v2_, dent_, Î¦2_] = {tableNroue[[1]], tableNroue[[2]],
tableNroue[[3]]};

In[128]:=
F[u1_,v1_,u2_,v2_,Î¦2_,k_]={sp11[u1,v1,dent,Î¦1][[1]]-sr22[u2,v2,dent,Î¦2][[1]],
sp11[u1,v1,dent,Î¦1][[2]]-sr22[u2,v2,dent,Î¦2][[2]],
sp11[u1,v1,dent,Î¦1][[3]]-sr22[u2,v2,dent,Î¦2][[3]],
np11[u1,v1,dent,Î¦1][[1]]-k*nr22[u2,v2,dent,Î¦2][[1]],
np11[u1,v1,dent,Î¦1][[2]]-k*nr22[u2,v2,dent,Î¦2][[2]],
np11[u1,v1,dent,Î¦1][[3]]-k*nr22[u2,v2,dent,Î¦2][[3]]};

In[129]:=
var={u1,v1,u2,v2,Î¦2,k};

(*here is the jacobian matrix*)

tic1 = AbsoluteTime[];
For[j = 1, j â?¤ Length[var], j++, gradtempo[[1, j]] = D[F[u1, v1, u2,
v2, Î¦2,
k][[1]], var[[j]]]];
For[j = 1, j â?¤ Length[var],
j++, gradtempo[[2, j]] = D[F[u1, v1, u2, v2, Î¦2, k][[2]],
var[[j]]]];
For[j = 1, j â?¤ Length[var], j++, gradtempo[[3, j]] = D[F[u1, v1, u2,
v2, Î¦2,
k][[3]], var[[j]]]];
For[j = 1, j â?¤ Length[var],
j++, gradtempo[[4, j]] = D[F[u1, v1, u2, v2, Î¦2, k][[4]],
var[[j]]]];
For[j = 1, j â?¤ Length[var], j++, gradtempo[[5, j]] = D[F[u1, v1, u2,
v2, Î¦2,
k][[5]], var[[j]]]];
For[j = 1, j â?¤ Length[var],
j++, gradtempo[[6, j]] = D[F[u1, v1, u2, v2, Î¦2, k][[6]],
var[[j]]]];

grad[u1_, v1_, u2_, v2_, Î¦2_, k_] = gradtempo;
g[u1_, v1_, u2_, v2_, Î¦2_, k_] = Norm[F[u1, v1, u2, v2, Î¦2, k]]^2;
tac1 = AbsoluteTime[];
tac1 - tic1

(*and here is the inverse of this matrix, that i would like to be
faster*)

tic = AbsoluteTime[];
Inverse[grad[0.5, 0.5, 0.5, 0.5, 0, 1];];
tac = AbsoluteTime[];
tac - tic

Sorry for my bad english

Best Regards,
JP

```

• Prev by Date: Intercepting and Controlling Mathematica Exceptions.
• Next by Date: Re: SubValues
• Previous by thread: Intercepting and Controlling Mathematica Exceptions.
• Next by thread: Re: SubValues