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]:=
(*Définition de la géométrie*)
(*Définition des géométries nominales dans R7 et R8*)
(*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*)
Ω[v1_]=27*omegamaxi/4*(v1-Rp)^2/(Rg-Rp)^2*(1-(v1-Rp)/(Rg-Rp))//Expand;
(*Géométrie du pignon avec bombé dans R7*)
sb1[u1_,v1_]={x1[u1,v1],y1[u1,v1]*Cos[Ω[v1]]+
z1[u1,v1]*Sin[Ω[v1]],-y1[u1,v1]*Sin[Ω[v1]]+z1[u1,v1]*Cos[Ω[v1]]}//\
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;
(*Calcul des normales des surfaces déformées*)
(*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 \
résolution*)
(*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]:=
(*Résolution*)
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};
gradtempo=Table[i-j,{i,1,Length[F[u1,v1,u2,v2,Φ2,k]]},{j,1,Length[var]}];
(*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