<<PhysicalConstants`
g=AccelerationDueToGravity/.{Meter->1,Second->1}//N;

Matematicas

Cálculo variacional

Problema libre
W==(∫_(Subscript[p, 0])^(Subscript[P, 1]))L(x,y(x),y'(x))dx
((∂L)/(∂y))-(∂/(∂x))((∂L)/(∂y'))==0

Identidad Beltrami 
L-y'((∂L)/(∂y'))==C

Problema constreñido
W=(∫_(Subscript[P, 0])^(Subscript[P, 1]))L(x,y(x),y'(x))dx, (∫_(Subscript[P, 0])^(Subscript[P, 1]))y(x)dx=C

1. Ejemplo general:  L(x,y(x),y'(x)):=y'(((x))^2)+3 y(x)x, y(0)=0, y((π/2))=1

Clear[x,f,df,y]
L[x_,f_,df_]:=df^2+3 x f
EDE[L_]:=D[L,y[x]]-D[D[L,y'[x]],x];
eq=EDE[L[x,y[x],y'[x]]]
s=DSolve[{eq==0,y[0]==0.0,y[π/2]==1},y[x],x]//Flatten
ff=y[x]/.s;
dff=D[y[x]/.s,x];
Integrate[L[x,ff,dff],{x,0.0,π/2}]
g1=Plot[ff,{x,0,Pi/2},PlotStyle->Thick];
g2=ListPlot[{{{0,0},{π/2,1}},{{0,0},{π/2,1}}},PlotStyle->{{Red,PointSize[0.02]},{White,PointSize[0.015]}}];
Show[g1,g2]

2. Braquistocrona curva

Clear[x,f,df,y,g,k,ξ,η,θ];
L[x_,f_,df_]:=Sqrt[1+df^2]/Sqrt[2g f];
BDE[L_]:=L-y'[x]D[L,y'[x]]
eq=BDE[L[x,y[x],y'[x]]]==k//FullSimplify
ξ[θ_]:=1/2k^2(θ-Sin[θ]);
η[θ_]:=1/2 k^2(1-Cos[θ]);
(1+(D[η[θ],θ]/D[ξ[θ],θ])^2)*η[θ]//FullSimplify
ParametricPlot[{ξ[θ],-η[θ]}/.k->1,{θ,0,3π/4},PlotStyle->{Red,Thick}]

3. Catenaria

Clear[x,f,df,ρ, g,u]
L[x_,f_,df_]:=ρ g f*Sqrt[1+df^2]+λ Sqrt[1+df^2];
eq=BDE[L[x,y[x],y'[x]]]//FullSimplify
sol=DSolve[eq==k,y[x],x]
Plot[y[x]/.sol[[1]]/.λ->0/.g->1/.ρ->1/.C[1]->1/.k->2,{x,-2,2},PlotStyle->Thick]

4.  Superfície mínima

Clear[x,y,f,df,ρ, g,u]
L[x_,f_,df_]:=2π f*Sqrt[1+df^2];
eq=BDE[L[x,y[x],y'[x]]]//FullSimplify
sol=DSolve[{eq==k},y[x],x]/.C[1]->0/.k->1//Flatten
ParametricPlot3D[{x,y[x]*Cos[η],y[x]*Sin[η]}/.sol[[1]],{x,0,0.3},{η,0,2Pi}]

Física

El movimiento está descrito por la física de Galileo. Se compone de una sola sentencia: todo movimiento minimiza el cambio. En la naturaleza, el cambio se mide por la acción física W.
Más precisamente, el cambio se mide por la diferencia de tiempo promedio entre la energía cinética T y energía potencial V. En otras palabras, el movimiento obedece al asi llamado principio de acción minima:

	δ W==0 , donde W==∫(T-V)dt

	L==T-V==L(x,(OverDot[x]),t),   δ W== ((∂L)/(∂x))-(∂/(∂t))((∂L)/(∂(OverDot[x])))==0

En otros términos, la naturaleza es tan floja como sea posible. La pereza natural del movimiento cotidiano sigue siendo válida en toda la física moderna, agregando unas pocas sentencias de 
límite, como v<c (donde c es la velocidad de luz) en la teoria de relatividad especial y W≥ℏ 
donde (ℏ es la constante de Planck reducida) en la mecanica cuantica etc...

Acción mínima

1. Tiro inclinado

Energia potencial V=g y(t)
Energia cinética T=(1/2) m (x'(((t))^2)+(y'(((t))^2))
Sistema de ecuaciones de Euler  ((∂L)/(∂x(t)))-(∂/(∂t))((∂L)/(∂x'(t)))==0,((∂L)/(∂y(t)))-(∂/(∂t))((∂L)/(∂y'(t)))==0

Clear[x,y,t]
V=g y[t];
T=m/2(D[x[t],t]^2+D[y[t],t]^2);
L=T-V;
δW={D[L,x[t]]-D[D[L,D[x[t],t]],t]==0,
D[L,y[t]]-D[D[L,D[y[t],t]],t]==0};
sol=Flatten[DSolve[δW,{x[t],y[t]},t]]

Condiciones iniciales

Off[Solve::ratnz]
Clear[k]; (* k cantidad de polvora *)
h=10;(* posicion de meta *)
m=1; (* masa de proyectilo *)
cond={x[0]==0,y[0]==0,x'[0]==k,y'[0]==k}; 
(* condiciones iniciales, angulo de disparo 45 grados*)
sol1=DSolve[Join[δW,cond],{x[t],y[t]},t]//Flatten (* solver el sistema de ode *)
sol2=Solve[{(y[t]/.sol1)==0,(x[t]/.sol1)==h},{t,k}]//Last (* tiro en blanco *)
g1=ParametricPlot[{x[t],y[t]}/.sol1/.t->τ/.sol2,{τ,0,0.1+t/.sol2}];
p={{x[t],y[t]}}/.sol1/.sol2;
g2=ListPlot[{p,p},PlotStyle->{{Red,PointSize[0.02]},{White,PointSize[0.015]}}];
Show[g1,g2]

Fricción

T=2;
μ=0.4; (* coeficiente de friccion *)
kk=9;  (* *)
δW={D[L,x[t]]-D[D[L,D[x[t],t]],t]==μ Abs[D[x[t],t]]D[x[t],t],
D[L,y[t]]-D[D[L,D[y[t],t]],t]==μ Abs[D[y[t],t]]D[y[t],t]};
(* sestema de ode *)
sol1=NDSolve[Join[δW,{x[0]==0,y[0]==0,x'[0]==k,y'[0]==k}/.k->kk],{x[t],y[t]},{t,0,T}]//Flatten; (* resolver el sistema *)
at=FindRoot[y[t]/.sol1,{t,T}];
a=x[t]/.sol1/.at
g1=ParametricPlot[{x[t],y[t]}/.sol1,{t,0,T}];
g2=ListPlot[{{{a,0}},{{a,0}}},PlotStyle->{{Red,PointSize[0.02]},{White,PointSize[0.015]}}];
Show[g1,g2]

Dar en blanco

d=6;
Clear[u]
alcance[k_]:=Block[{x,y,t,μ=0.4,T=2,s1,s2},
s1=NDSolve[{-x''[t]==μ Abs[x'[t]]x'[t],-g-y''[t]==μ Abs[y'[t]]y'[t],x[0]==0,y[0]==0,Derivative[1][x][0]==k,Derivative[1][y][0]==k},{x[t],y[t]},{t,0,T}]//Flatten;
s2=FindRoot[y[t]==0/.s1,{t,T}];
Return[x[t]/.s1/.s2]]

FindRoot[alcance[u]==d,{u,10},Evaluated->False]

2. Doble péndulo

La función de Lagrange

s={m1->1,m2->1,l1->1,l2->1,g->9.8};
x1[t_]:= l1 Sin[θ1[t]];
y1[t_]:=- l1 Cos[θ1[t]];
x2[t_]:=(l1 Sin[θ1[t]]+ l2 Sin[θ2[t]]);
y2[t_]:=-(l1 Cos[θ1[t]]+ l2 Cos[θ2[t]]);

V=m1 g y1[t]+m2 g y2[t];
T=1/2 m1 (D[x1[t],t]^2+D[y1[t],t]^2)+1/2 m2 (D[x2[t],t]^2+D[y2[t],t]^2)//Expand//FullSimplify;
L=(T-V)//FullSimplify;
eq={D[L,θ1[t]]-D[D[L,D[θ1[t],t]],t]==0,
D[L,θ2[t]]-D[D[L,D[θ2[t],t]],t]==0,θ1[0]==0,θ2[0]==2,θ1'[0]==0,θ2'[0]==0}/.s//Expand//FullSimplify;

Solución de la equación diferencial

f=NDSolve[eq,{θ1[t],θ2[t]},{t,0,6}];

Gráfica

g1=ParametricPlot[{x1[t],y1[t]}/.f/.s,{t,0,5}];
g2=ParametricPlot[{x2[t],y2[t]}/.f/.s,{t,0,5}];

h0=Table[{0,0},{t,0,5,0.1}];
h1=Flatten[Table[{x1[t],y1[t]}/.f/.s,{t,0,5,0.1}],1];
h2=Flatten[Table[{x2[t],y2[t]}/.f/.s,{t,0,5,0.1}],1];
g3=Map[Graphics[Line[#],PlotRange->{{-2,2},{-2,0}}]&,Transpose[{h0,h1,h2}]];
Show[g1,g2,PlotRange->All]

Animación

ListAnimate[g3]

3. Oscilador harmónico clásico

Lagrange : L(t,x,x')==(1/2) m ((x')^2)-((k (x^2))/2)
Hamilton : H(t,q,p)=((k (q^2))/2)+((p^2)/(2 m))
Ecuaciones : ((∂H)/(∂p))==((∂q)/(∂t)), ((∂H)/(∂q))==-((∂p)/(∂t))

Clear[p,q,t,m,k]
T=1/2 p[t]^2/m;
V=1/2 k q[t]^2;
H=T+V

eq={D[H,p[t]]==D[q[t],t],D[H,q[t]]==-D[p[t],t]}
sol1=FullSimplify[(DSolve[Join[eq,{q[0]==0,p[0]==1}],{q[t],p[t]},t]/.k->ω^2 m/.m->1),Assumptions->{m>0,k>0}]

Probabilidad

Clear[x,P,A]
x[t_]=A Cos[ω t];
τ[ξ_]=ArcCos[ξ/A]/ω;
v[ξ_]=-D[x[t],t]/.t->τ[ξ];
P[ξ_]=ω/Pi/v[ξ]
Integrate[P[ξ],{ξ,-A,A},Assumptions->A>0];
Plot[P[ξ]/.A->1,{ξ,-1,1},PlotStyle->Thickness[0.006],PlotRange->{0,1}]

4. Péndulo matemático

Hamiltonián  H=g l m (1-cos(ϕ(t)))+(1/2) m (((l ψ(t)))^2)
ϕ(t), ψ(t)=ϕ'(t)

Clear[H,ϕ,ψ,t,m,l]
H=m*(l*ψ[t])^2/2+m g l (1-Cos[ϕ[t]]);
eqn={D[H,ψ[t]]==D[ϕ[t],t],D[H,ϕ[t]]==-D[ψ[t],t]}/.l->1/.m->1;
ψ0=6.26;
T0=15;
sol=NDSolve[Join[eqn,{ϕ[0]==0,ψ[0]==ψ0}],{ϕ[t],ψ[t]},{t,0,T0}];
Plot[ϕ[t]/.sol,{t,0,T0}]
ParametricPlot[{ϕ[t],ψ[t]}/.sol,{t,0,T0}]

Diagrama de fase

ψ0=6.3;
T0=5;
sol=NDSolve[Join[eqn,{ϕ[0]==0,ψ[0]==ψ0}],{ϕ[t],ψ[t]},{t,0,T0}];
Plot[ϕ[t]/.sol,{t,0,T0}]
ParametricPlot[{ϕ[t],ψ[t]}/.sol,{t,0,T0}]

5. Péndulo elástico

Hamiltonián  H==(1/(2m))p.p+m q.g+(k/2) (((L-(Sqrt[q.q])))^2)

m=1; k=100;L=1; T=5;
Clear[px,py,qx,qy,p,q,t]
VD[U_,Q_]:=Map[D[U,#]&,Q]
p[t_]={px[t],py[t]};
q[t_]={qx[t],qy[t]};
H=p[t].p[t]/2m+m q[t].{0,g} + k (L-Sqrt[q[t].q[t]])^2/2
eqns={MapThread[#1==#2&,{VD[H,p[t]],D[q[t],t]}],MapThread[#1==#2&,{VD[H,q[t]],-D[p[t],t]}]}//Flatten//FullSimplify;
inicond={px[0]==1,py[0]==0,qx[0]==0,qy[0]==-L};
sol=NDSolve[Join[eqns,inicond],Join[p[t],q[t]],{t,0,T}];
g1=ParametricPlot[{qx[t],qy[t]}/.sol,{t,0,T}]

Animación

Animate[Show[g1,ListPlot[{{qx[t],qy[t]}}/.sol/.t->τ,PlotRange->{{-0.4,0.4},{-1.3,-0.90}},PlotStyle->{Red,PointSize[0.02]}]],{τ,0,T}]

Ecuaciones diferenciales parciales

1. Oscilaciones radiales propias de la membrana circular

(((∂^2)U(r,t))/(∂(r^2)))+(1/r)((∂U(r,t))/(∂r))==(1/(c^2))(((∂^2)U(r,t))/(∂(t^2))), |U(0,t)|<∞, |U(a,t)|=0
U(r,t)==R(r) T(t)
((R'(t)+r R''(t))/(r R(t)))=((T(t))/((c^2) T(t)))=-(k^2)

Clear[U,R,T,c,k,n,A,B,s,a]
U[r_,t_]:=R[r]*T[t]
eq=Map[#==-k^2&,({D[U[r,t],{r,2}]+1/r D[U[r,t],r],1/c^2*D[U[r,t],{t,2}]}/U[r,t]//FullSimplify)];

de={R[r],T[t]};
nd={r,t};
sol=(MapThread[DSolve[#1,#2,#3]&,{eq,de,nd}]//Flatten)

sx=sol[[1]]/.C[2]->0/.C[1]->1;
sy=sol[[2]]/.C[1]->A[n]/.C[2]->B[n];
s[n_]=R[r]*T[t]/.sx/.sy/.k->BesselJZero[0,n]/a

Gráfica

Clear[f]
a=1;
c=1;
f[r_,t_]=(s[4]/.A[4]->1/.B[4]->0)//N
ξ=BesselJZero[0,1]//N
Animate[ParametricPlot3D[{r Cos[φ],r Sin[φ],f[r,t]},{r,0,a},{φ,0,2 π},PlotPoints->30,PlotRange->{{-a,a},{-a,a},{-1,1}}],{t,0,(2 π a)/(c ξ)}]

2. Osclaciones del hilo colgado

μ (((∂^2)u(x,t))/(∂(t^2)))=(∂/(∂x))(τ(x) ((∂u(x,t))/(∂x))), 
|u(0,t)|<∞, |u(l,t)|==0, u(x,0)=L-x
Tensión de hilo: τ(x)==g μ x ⇒ (c^2)==g 
(1/(c^2))(Subscript[u, tt])=(Subscript[((Subscript[xu, xx])), x]),⇒ x=(ξ^2)  ⇒(2/(c^2))(Subscript[u, tt])=(1/ξ)(Subscript[((Subscript[ξu, ξ])), ξ])
u(ξ,t)=Ξ(ξ)T(t)⇒ Ξ''+(1/ξ)Ξ'λ Ξ(ξ)==0, T''+((((2/c)))^2)λ T=0
(Subscript[u, n])(x,t)=(Subscript[J, 0])((Subscript[ξ, n])(Sqrt[(x/l)]))((Subscript[A, n]) cos(t (Subscript[ω, n]))+(Subscript[B, n]) sin(t (Subscript[ω, n])))

Clear[x, t, ξ, λ, ω, c, L, u]
L = 1; c = 1;
f[x_] := (L - x);
M = 8;
λ = 
    Table[BesselJZero[0, n]/Sqrt[L], {n, 1, M}] // N;
ω = (1/2) c λ;
A = 
    Integrate[ξ f[ξ^2] BesselJ[0, λ ξ], {ξ, 0, 
          Sqrt[L]}]/
      Integrate[ξ  BesselJ[0, λ ξ]^2, {ξ, 0, 
          Sqrt[L]}];
u[ξ_, t_] = (A Cos[ω t]).BesselJ[0, λ ξ];
Animate[Rotate[Plot[u[Sqrt[x],t],{x,0,L},PlotStyle->Thick,PlotRange->{{0,L},{-1.2,1.2}},Axes->False],Pi/2],{t,0,10,0.001}]

3. Oscilaciones propias de la cuerda incrustada

(((∂^2)u(x,t))/(∂(t^2)))==(c^2) (((∂^2)U(x,t))/(∂(x^2))), u(0,t)==0, u(L,t)==0
(Subscript[u, n])(x,t)==(Subscript[A, n]) cos(t (Subscript[ω, n])) sin((Subscript[λ, n]) x), (Subscript[λ, n])==((π n)/L), (Subscript[ω, n])=c (Subscript[λ, n])

L=1;c=1;n=3;
λ=n π/L;
ω=c λ;
u[x_,t_]:=Cos[ω t] Sin[λ x];
Animate[Plot[u[x,t],{x,0,L},PlotStyle->Thick,PlotRange->{{0,L},{-1,1}}],{t,0,(2 L)/c,(2 L)/(100 c)}]

4. Oscilaciones generales de la cuerda incrustada

u(x,0)=f(x), (Subscript[u, t])(x,0)==0

Clear[n,λ, ω,t,x,u,f,g]
L=π;c=1; M=10;
λ[n_]:=n π/L;
ω[n_]:=c λ[n];
u[x_,t_,n_]:=Cos[ω[n] t] Sin[λ [n]x]
g[x_]:=Table[u[x,0,n],{n,M}]
f[x_]:=Piecewise[{{1,0.2<x<0.6}},0]
A=2/L Integrate[f[x]*g[x],{x,0,L}]
g1=Plot[f[x],{x,0,π},PlotStyle->{Thick,Red}];
g2=Plot[Evaluate[A.g[x]],{x,0,L}];
Show[g1,g2,PlotRange->All]

Animación

Animate[Plot[Evaluate[A.Table[u[x,t,n],{n,M}]],{x,0,L},PlotStyle->Thick,PlotRange->{{0,L},{-1,1}}],{t,0,(2 L)/c,(2 L)/(100 c)}]

5. Cuántum hármonic oscilator

(1/(2m))(((∂^2)ψ(x))/(∂(x^2)))+V (x)ψ(x)=E ψ(x), V(x)=((m (ω^2)(x^2))/2)
(Subscript[ψ, n])(x)=(e^(-(x^2)/2))(Subscript[H, n])(x)/(2^nn!(Sqrt[π]))

Block[{n=20},Plot[{(E^-x^2 HermiteH[n,x]^2)/(2^n n! Sqrt[π]),Piecewise[{{1/(Sqrt[2 n+1-x^2] π),x^2<2 n+1}}]},{x,-Sqrt[3 n],Sqrt[3 n]},PlotPoints->200,PlotRange->{0,0.3}]]

Relatividad

1. Invariante

Clear[U,v,c,u,w]

b={β->v/c,γ->1/Sqrt[1-β^2]};
L={{γ,-γ β},{-γ β,γ}};

u[t_,x_]=(L.{c t,x})/{c,1}//.b//FullSimplify
LI=Inverse[L]
w[τ_,ξ_]=LI.{c τ,ξ}/{c,1}//.b//FullSimplify

D[U[{τ,ξ}],{τ,2}]-c^2 D[U[{c τ,ξ}],{ξ,2}]//FullSimplify
D[U[w[τ,ξ]],{τ,2}]-c^2 D[U[w[τ,ξ]],{ξ,2}]//FullSimplify