Ejemplos para Física 1

Oscilador anarmonico

Planteamos las ecuaciones de Newton en general

Eqs:=m x''[t] == F[x[t]] (* Ecuacion de Newton en general para una dimension *)  
F[x_]:= -a (x^3)+b x (* Fuerza del oscilador anarmonico *)

Si la fuerza es la elastica del resorte (a=0 y b<0)

a=0;b=-1;m=1; (* Elijo algunos valores de los parametros *)
sol=DSolve[{Eqs,x[0]==x0,x'[0]==v0},x[t],t] //Flatten(* Resuelve la ecuacion diferencial analiticamente, si puede... *)

Ecuaciones Horarias

Plot[Evaluate[{x[t]/.sol,D[x[t]/.sol,t],D[x[t]/.sol,{t,2}]}/.{x0->1,v0->0}],{t,0,20},AxesLabel-> {"t","x,v"}]

Espacio de Fases

ParametricPlot[Evaluate[{x[t]/.sol,D[x[t]/.sol,t]}/.{x0->1,v0->0}],{t,0,50},AspectRatio->1,AxesLabel-> {"x","v"}]

Caso Anarmonico

a=10;b=-1;m=1; (* Elijo nuevos valores, ahora con a ≠ 0 para que sea anarmonico *) 
sola=NDSolve[{Eqs,x[0]==1,x'[0]==0},x[t],{t,0,100}]//Flatten

Ecuaciones Horarias

Plot[Evaluate[{x[t]/.sola,D[x[t]/.sola,t]}],{t,0,20},AxesLabel-> {"t","x,v"}]

Espacio de Fases

ParametricPlot[Evaluate[{x[t]/.sola,D[x[t]/.sola,t]}],{t,0,20},AspectRatio->1,AxesLabel->{"x","v"}]

a=1;b=1;m=1; (*otro  ejemplo, cualitativamente distinto, por que? *)
sola=NDSolve[{Eqs,x[0]==1.6,x'[0]==0},x[t],{t,0,100}]//Flatten;
Plot[Evaluate[{x[t]/.sola,D[x[t]/.sola,t]}],{t,0,60},AxesLabel-> {"t","x,v"},AspectRatio->0.3]

Energia (como en el ejercicio de la guia)

DSolve[{-V'[x]== F[x],V[0]==0},V[x],x] (* Utiliza los valores en memoria de a y b, si ya estan definidos *)

V[x_]:= a ((x^4)/4)- b ((x^2)/2) (* Potencial en general, el := hace que se actualice automaticamente si cambia a o b *)
Em[x_,v_]:=(1/2)m (v^2)+ V[x] (* Energia Mecanica total en funcion de (x,v) *)
Xeq[a_,b_]:= Solve[D[V[x],{x,1}]==0,x]; (* Puntos de equilibrio, podria escribirlo directamente en terminos de F[x], claro *)

Tratamiento Analitico

a=.;b=.;xeq=Xeq[a,b]

seri=Table[Series[V[x],{x,x/.xeq[[i]],2}],{i,1,3}] (* Calculo el desarrollo en Serie hasta orden 2 alrededor de los puntos de equilibrio *)

Coefficient[seri,x,2] (* De cada elemento de la lista tomo el coeficiente que multiplica (x^2) , que son k/2 del oscilador aproximado de peq. oscilaciones *)

(* Diagramas de Potencial *)
a=1;b=1;plot1=Plot[V[x],{x,-2,2},AxesLabel->{"x","V[x]"},Epilog-> {Text["a="<> ToString[a],{0.4,0.5}],{RGBColor[1,0,0],Line[{{-2,#},{2,#}}]& /@ {1,2,0,-0.2}}}]

Xeq[a,b]

V[x]/.Xeq[a,b]

plot1b=ContourPlot[Evaluate[{Em[x,v]==1, Em[x,v]==2,Em[x,v]==0,Em[x,v]==-0.2}],{x,-3,3},{v,-3,3}]

Como hacemos para calcular el período?
De la conservación de la energía tenemos  que E=(1/2)(mv^2)+V[x]  no cambia, entonces la velocidad como función de la posición queda v[x] = (Sqrt[(2 (E - V[x]))/m])   y por lo tanto como v[x]=(dx/dt)   integrándola entre dos puntos de retorno me da la mitad del período,  (Overscript[(∫_(pto de (Subscript[retorno, 1]))), pto de (Subscript[retorno, 2])])(dx/(v[x]))=(Overscript[(∫_0), τ/2])dt.  En notación del Mathematica sería:

velocidad[x_,Em_]:=(Sqrt[((2(Em- V[x]))/m)]) (* Modulo de la velocidad en una dada posicion*)
Periodo[Em_,pts__:{1,2}]:=Block[{Em0=Em},
retorno=Select[(x/. Solve[velocidad[x,Em0]==0,x]),# ∈ Reals &];NIntegrate[(1/(velocidad[x,Em0])),{x,retorno[[pts[[1]]]],retorno[[pts[[2]]]]}]2]

period1=Table[{Em,Periodo[Em]//Chop},{Em,-0.24999,1,0.01}];

ListPlot[Re[period1],Joined->True,AxesLabel->{"E","τ"},Epilog->{{RGBColor[1,0,0],Line[{{-0.25,((2π)/(Sqrt[2b]))},{1,((2π)/(Sqrt[2b]))}}]}}]

Como queda si la expresamos en terminos de la frecuencia y la amplitud inicial???.

a=1;b=-1;  (* Tomo ahora el mismo valor del principio de la practica *) plot2=Plot[V[x],{x,-2,2},AxesLabel->{"x","V[x]"},Epilog->{Text["a="<>ToString[a],{0.4,0.5}],{RGBColor[1,0,0],Line[{{-2,#},{2,#}}]&/@{1,2}}}]

Xeq[a,b]

V[x]/.Xeq[a,b]

plot2b=ContourPlot[Evaluate[(Em[x,v]==#1&)/@{1,2}],{x,-3,3},{v,-3,3}]

period2=Table[{Em,Periodo[Em]//Chop},{Em,0.0001,40,0.5}];

ListPlot[Re[period2],Joined->True,PlotRange->All,AxesLabel->{"E","τ"}]

a=-1;b=1;plot3=Plot[V[x],{x,-2,2},AxesLabel->{"x","V[x]"},Epilog->{Text["a="<>ToString[a],{0.4,0.5}],{RGBColor[1,0,0],Line[{{-2,#},{2,#}}]& /@ {1,0,-0.2}}},PlotRange->{-6,2}]

Xeq[a,b]

V[x]/.Xeq[a,b]

plot3b=ContourPlot[Evaluate[(Em[x,v]==#1&)/@{1,0,-0.2}],{x,-3,3},{v,-3,3}]

a=-1;b=-1;plot4=Plot[V[x],{x,-2,2},AxesLabel->{"x","V[x]"},Epilog->{Text["a="<>ToString[a],{0.4,0.5}],{RGBColor[1,0,0],Line[{{-2,#},{2,#}}]& /@ {-0.2,0,0.2}}}]

Xeq[a,b]

V[x]/.Xeq[a,b]

plot4b=ContourPlot[Evaluate[(Em[x,v]==#1&)/@{0.2`,0,-0.2`}],{x,-3,3},{v,-3,3}]

period4=Table[{Em,Periodo[Em,{2,3}]//Chop},{Em,0.0001,(1/4),0.001}];

ListPlot[Re[period4],Joined->True]

RESUMEN DEL POTENCIAL ANARMONICO

Show[GraphicsGrid[{{plot1,plot2,plot1b,plot2b},{plot3,plot4,plot3b,plot4b}}]]

Gravitacion - Problema de Kepler (para masas iguales)

Como vieron en teorica el problema tiene solucion analitica sencilla separando el problema en el centro de masa y movimiento relativo. Aca lo que veremos es como resolverlo numericamente (obviamente).

Ecuaciones de Newton en dos dimensiones y coordenadas cartesianas.

(Subscript[F, 1])=m (Subscript[a, 1  ]) y  (Subscript[F, 2]) = m (Subscript[a, 2]) con  (Subscript[F, 1]) y (Subscript[F, 2]) las fuerzas gravitatorias que siente cada uno. 
(Subscript[F, 1])= - ((G m m ((Subscript[r, 1])-(Subscript[r, 2])) )/(|(Subscript[r, 1])-(Subscript[r, 2])((|)^3)))   y  (Subscript[F, 2])= - ((G m m ((Subscript[r, 2])-(Subscript[r, 1])) )/(|(Subscript[r, 2])-(Subscript[r, 1])((|)^3)))
Estas ecuaciones diferenciales en 2D  en notacion del Mathematica se escriben

eqG:={ x1''[t]== -Gm (((x1[t]-x2[t]))/((((Sqrt[(((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)])))^3)),y1''[t]==-Gm (((y1[t]-y2[t]))/((((Sqrt[(((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)])))^3)),x2''[t]==-Gm (((x2[t]-x1[t]))/((((Sqrt[(((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)])))^3)),y2''[t]==-Gm (((y2[t]-y1[t]))/((((Sqrt[(((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)])))^3)),x1[0]==-d/2,x2[0]==d/2,y1[0]==0,y2[0]==0,x1'[0]==v0 Cos[α], y1'[0]== v0 Sin[α],x2'[0]== -v0 Cos[α],y2'[0]==-v0 Sin[α]};

Ya puse como condicion inicial  una relacion entre las velocidades iniciales (Problema 7. Guía de Gravitación )

E0:=(v0^2) -(Gm/d) (* Energia inicial  / m *)

(Subscript[L, cm])=m (Subscript[r, 1])×(Subscript[v, 1])+m (Subscript[r, 2]) × (Subscript[v, 2])= μ r × vdonde μ es la masa reducida, y  r  y v  son las coordenadasa relativas.

Lcm:= v0 d Sin[α]  (* Momento angular total inicial / m *)

μ=1/2 ;(* Masa reducida (por unidad de masa) *)

Energia + Momento angular  (* todo dividido por m*)

Em[r_,vr_]:= (1/2)μ (vr^2)+Veff[r]
Veff[r_]:=(1/2)((Lcm^2)/(μ (r^2)))-((Gm )/r)

Resolvemos

Gm=1;d=0.75;v0=1;α=30*Pi/180;tmax=100; (* Elijo condiciones iniciales y otros parametros *) sol=NDSolve[eqG,{x1[t],y1[t],x2[t],y2[t]},{t,0,tmax}]//Flatten;
vx1[t]=D[x1[t]/.sol,t];vy1[t]=D[y1[t]/.sol,t];vx2[t]=D[x2[t]/.sol,t];vy2[t]=D[y2[t]/.sol,t];
(* Grafico las trayectorias *)
ParametricPlot[Evaluate[{{x1[t],y1[t]},{x2[t],y2[t]}}/.sol],{t,0,tmax},Frame->True,PlotStyle->{{Thickness[0.003],RGBColor[1,0,0]},{Thickness[0.003],RGBColor[0,1,0]}},AspectRatio->1,AxesLabel->{"x","y"},Prolog->{PointSize[Large],Point[{{x1[t],y1[t]},{x2[t],y2[t]}}/.sol/.t->0]}]

Plot[Veff[r],{r,0,3},Epilog->{{RGBColor[1,0,0],Line[{{0,E0},{3,E0}}]},Text["(Subscript[L, cm])="<> ToString[Lcm],{1,-0.5}]}]

Movimiento relativo

ParametricPlot[{x1[t]-x2[t],y1[t]-y2[t]}/.sol,{t,0,15},Frame->True,PlotStyle->{Thickness[0.003],RGBColor[1,0,0]},AspectRatio-> 1,AxesLabel-> {"x","y"}]

Leyes de Kepler

(τ^2)=4 ((π^2)/Gm) (r^3) (* τ es el período y r es el semieje mayor *)

Como lo encontrarias aca ??

Magnitudes conservadas

Impulso Angular

L[x_,v_]:=m x \[Cross]v;
L1=L[{x1[t],y1[t],0},{vx1[t],vy1[t],0}].{0,0,1}/.sol; (* Componente z del momento angular de 1, todo *)
L2=L[{x2[t],y2[t],0},{vx2[t],vy2[t],0}].{0,0,1}/.sol;(* Componente z del momento angular de 2, todo *)

Plot[L1,{t,0,tmax}] (* Deberia ser constante, pero.... *)

Energia Mecánica Total

Ec1=(1/2)m (((vx1[t])^2)+((vy1[t])^2));Ec2=(1/2)m(((vx2[t])^2)+((vy2[t])^2));Epot12=-(Gm/(Sqrt[(((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)]))/.sol;

Plot[{Ec1,Ec2,Epot12,Ec1+Ec2+Epot12},{t,0,tmax},AxesLabel->{"t","Energias"}]