(* TUTORIAL  (Para Ejecutar instrucción: shift+enter o Intro en num. pad) *)

1+1

(* guardar resultado en a *)

a=1+1

(* al finalizar con ; se ejecuta la instruccion pero no se muestra el resultado (importante para salidas grandes) *)

a=1+1;

(* Suma, resta, multiplicacion, division, raiz, etc.*)

{2+3,2-3,2*3,2  3, 2/3,2./3,2^3,Sqrt[2],Sqrt[2.],2^(1/2), 2 ^(1./2)}

(* Sqrt[2] con k decimales: N[Sqrt[2],k] *)

N[Sqrt[2],1000]

{Pi,E}

N[Pi,1000]

N[E,1000]

N[EulerGamma,100]

100!

N[100!]

N[100!,50]

(* es posible manejar numeros muy grandes: *)

N[1000!]

1000!

(* Ploteo de funciones, derivadas, etc. *)

Clear[x];

Plot[Sin[x],{x,-2 Pi,2 Pi},PlotStyle->Thick,AxesStyle->14,AxesLabel->{"x","sen(x)"}]

Plot[Sin[x],{x,-2 Pi,2 Pi},PlotStyle->Thick,AxesStyle->14,AxesLabel->{"x","sen(x)"},Filling->Axis]

Sin'[x]

D[Sin[x],{x,1}]

Exp[r*x]

D[Exp[r*x],x]

(* Derivada 22 *)

D[Exp[r*x],{x,22}]

Integrate[Sin[x],x]

∫Sin[x] \[DifferentialD] x

(* Integral definida *)

Integrate[Sin[x],{x,0,Sqrt[2]*Pi/2}]

Integrate[Exp[-x^2],x]

(* Integracion analitica Integrate; Integracion numerica: Nintegrate *)

NIntegrate[Exp[-x^2],{x,0,1}]

(* comprobacion *)

N[(1/2) (Sqrt[π]) Erf[1]]

Integrate[Exp[-x^p],x]

(* Formas equivalentes de evaluar funciones: f[x] o x//f *)

Sin[Pi/2]

Pi/2//Sin

Sin[x]

x//Sin

(* Definicion de funcion: Nombre[x_]:= Expresión de la función *)

f[x_]:=Sin[x^2/Tan[Sqrt[1+x^4/Exp[x]]]]

Plot[f[x],{x,-1,12},PlotStyle->Thick,AxesStyle->14,AxesLabel->{"x","f(x)"}]

f'[x]

Plot[{f[x],f'[x]},{x,-1,2},PlotStyle->Thick]

FindRoot[f'[x]==0,{x,1}]

x0=x/.%

x0

(* Serie de Taylor alrededor de a hasta orden n: Series[f[x],{x,a,n}]; 
para expresion normal sin O(x^(n+1)): N[Series[f[x],{x,a,n}] *)

Series[f[x],{x,0,2}]

Chop[N[%]]

p2[x_]=Series[f[x],{x,0,2}]

p4[x_]=Series[f[x],{x,0,2}]

p8[x_]=Series[f[x],{x,0,8}]

Chop[N[%]]

p2n[x_]=Normal[p2[x]];p8n[x_]=p8[x]//Normal

pl1=Plot[{f[x],p2n[x],p8n[x]},{x,-1,2},PlotStyle->{Thick,{Thick,DotDashed},{Thick,Dashed}},PlotRange->{{-1,2},{-1.1,1.1}},AxesStyle->14,AxesLabel->{"x","f(x)"}]

pl2=Show[pl1,Graphics[{Text[StyleForm["(Subscript[P, 2])(x)",16],{1.3,.75}],
Text[StyleForm["(Subscript[P, 8])(x)",{16}],{1.45,-.75}],Text[StyleForm["f(x)",{16}],{1.6,.4}]}]]

Series[Cos[x],{x,0,4}]

Normal[Series[Cos[x],{x,0,4}]]

Series[Cos[x],{x,0,20}]

Do[pol[n,x_]=Normal[Series[Cos[x],{x,0,n}]],{n,2,20,2}]

pol[2,x]

pol[4,x]

Clear[x,n];pl=Plot[{Cos[x],pol[2,x],pol[4,x],pol[20,x]},{x,-4*Pi,4*Pi},PlotStyle->{Thick,{Thick,Dashed},{Thick,Dotted},{Thick,Dashing[{0.02,0.01,0.01,0.01}]},
{Thickness[0.01],Dashed}},AxesStyle->14,AxesLabel->{"x","(Subscript[P, n])(x)"},PlotRange->{Automatic,{-2,2}}]

Clear[x,n];Show[pl,Graphics[{Text[StyleForm["(Subscript[P, 2])(x)",16],{4,-1.5}],
Text[StyleForm["(Subscript[P, 20])(x)",{16,Bold}],{11.5,1.75}]}]]

Series[Exp[-x^2],{x,0,14}]

(* Sumas: Sum[expresion,{n,valor inicial, valor final, paso}]  y Series *)

Sum[x^n,{n,1,10}]

Sum[x^n,{n,0,10,2}]

Sum[x^n,{n,0,Infinity,1}]

Sum[x^n,{n,0,Infinity,1},GenerateConditions->True]

Sum[x^n,{n,0,Infinity,2}]

Sum[x^n/n!,{n,0,Infinity}]

Sum[1/n^2,{n,1,Infinity}]

N[Sum[1/n^2,{n,1,Infinity}]]

a=%

N[Sum[1/n^2,{n,10}]]

N[Sum[1/n^2,{n,100}]]

N[Sum[1/n^2,{n,10000}]]

sm[m_]:=N[Sum[1/n^2,{n,m}]];

a\\\

ListPlot[{Table[{m,a},{m,100}],Table[{m,sm[m]},{m,100}]},PlotRange->{{0,100},{0.5,1.8}}]

Sum[1/n,{n,Infinity}]

Sum[1/n,{n,m}]

Series[%,{m,Infinity,2}]

smap[m_]=Normal[%]

N[EulerGamma,20]

N[smap[2]-HarmonicNumber[2]]

ListPlot[{Table[{m,HarmonicNumber[m]},{m,100}],
Table[{m,smap[m]},{m,100}]},PlotRange->{{0,100},{0.,10}}]

N[smap[10]-HarmonicNumber[10]]

smap[10.^23]

smap[10.^83]

Sum[(-1)^(n+1)/n,{n,Infinity}]

Clear[p];Sum[1/n^p,{n,Infinity}]

Zeta[2]

Plot[Zeta[p],{p,1,10},PlotRange->{{0,10},{0,5}},PlotStyle->Thick]

(* ejemplo de funcion definida mediante if *)

f=.;x=.;f[x_]:=If[x<0,Sin[x],2*Sin[x]+x/10];Plot[f[x],{x,-30,30},PlotStyle->Thick]

Plot3D[Sin[x*y],{x,-3,3},{y,-3,3}]

(* Vectores, Matrices, Autvalores (Eigenvalues) y Autovectores (Eigensystem) *)

(* para limpiar valores: Clear *)

Clear[a,b,c];

(* vectores: entre llaves; matrices: doble llave (vectores de vectores) *)

vec={a,b,c}

mat={{a,b},{c,d}}

(* para visualizar como matriz: MatrixForm *)
mat//MatrixForm

(* Producto escalar o matricial *)

vec.vec

mat.mat

mat//MatrixForm

(* Estos son productos elemento a elemento. No confundir *)

vec*vec

mat*mat

Tr[vec,Times]

mat={{1,1,1},{1,-1,2},{2,0,3}};

mat//mf

Transpose[mat]//MatrixForm

Det[mat]

Inverse[mat]

mat2={{1,1,1},{1,-1,2},{2,0,4}};

mat2inv=Inverse[mat2]

mat2.mat2inv

mat2.mat2inv//MatrixForm

(* Forma de Gauss Jordan (reduccion por filas)  *)

matr=RowReduce[mat]

matr//MatrixForm

(* fila 1 *)

mat[[1,All]]

(* Columna 1 *)

mat[[All,1]]

(* Rango *)

MatrixRank[mat]

(* Resolucion de sistema MX=B por 3 metodos equivalentes *)

Solve[{x+y+z==a,x-y+2z==b,2x+4z==c},{x,y,z}]

mat2//MatrixForm

Solve[mat2.{x,y,z}=={a,b,c},{x,y,z}]

Inverse[mat2].{a,b,c}

(* Resolucion de sistemas caso no singular  (CD) *)

mat//MatrixForm

matinv=Inverse[mat]

(* polinomio caracteristico,  Autovalores, autovectores *)

CharacteristicPolynomial[mat,x]

Det[mat-IdentityMatrix[3]*x]

(* Autovalores y autovectores *)

z=Eigenvalues[mat]

Solve[p==0,x]

v=Eigenvectors[mat]

(* autovalores y autovectores juntos: *)

Eigensystem[mat]

{z,v}=Eigensystem[mat]

z

v

v[[1]]

mat.v[[1]]

mat.v[[1]]-z[[1]] v[[1]]

s=Transpose[v]

d=Inverse[s].mat.s

d//MatrixForm

Inverse[s]

Inverse[s]//mf

Inverse[s].s//MatrixForm

(* Generacion de Matriz Random simétrica de nxn *)

n=4;ma=Table[Random[],{i,n},{j,n}];mas=(ma+Transpose[ma])/2

mas[[1]]

mas//MatrixForm

ma//MatrixForm

mf[au_]=MatrixForm[au];

mas//mf

Eigenvalues[ma]

Eigensystem[ma]

mas[[1,All]]

mas[[All,1]]

{z,v}=Eigensystem[ma];s=Transpose[v];si=Inverse[s];

z

s//MatrixForm

(* Determinante y traza *)

{Tr[a],Tr[z]}

{Det[a],Tr[z,Times]}

si=Inverse[s];s.si//MatrixForm//Chop

si.ma.s//MatrixForm//Chop

ma2=ma.ma;ma2//MatrixForm

ma2p=s.DiagonalMatrix[z^2].si;ma2p//mf

ma2-ma2p//MatrixForm//Chop

(* Gram Schmidt: partiendo de vec, gs genera vectores ortogonales en v y ortonormales (normalizados) en vn *)

Clear[vec,v];gs:=(v=vec;vn=v;le=Length[vec];Do[v[[i]]=vec[[i]]-If[i==1,0,Sum[(v[[j]].vec[[i]]/(v[[j]].v[[j]]))*v[[j]],{j,i-1}]];vn[[i]]=v[[i]]/Sqrt[v[[i]].v[[i]]],{i,le}]);

vec={{1,1,1},{1,1,2},{1,2,1}};

gs;v

v.Transpose[v]

vn

vn.Transpose[vn]

(* Orthogonalize: Instrucción que directamente genera vectores ortonormales *)

Orthogonalize[vec]

(* Ejercicios *)

"1) Encontrar los polinomios de Taylor alrededor del origen de la función 
 Exp[-x^2/10] Cos[x], de grado 2, 4,10, y 40;
Graficar la funcion y los poliniomios entre -2 Pi y 2 Pi;


2) Generar una matriz simétrica m de 6 x 6 con elementos aleatorios entre 
0 y 1  y diagonalizarla. Encontrar 
a) Sus autovalores,
b) Sus autovectores normalizados 
c) Una matriz s tal que st.m.s = d, con d diagonal, verificando el resultado. 
d) m^12, mediante dos métodos distintos, comprobando el resultado.  
e) Det[m], mediante dos métodos distintos, comprobando el resultado. 
f) Tr[m], mediante dos métodos distintos, comprobando el resultado. 
g) La inversa de m, en caso de que exista, mediante dos métodos distintos, comprobando el 
resultado. 

3) Idem 3) para una matriz m de 50 x 50 construida del igual forma. 

4) a) Generar 10 vectores de (Superscript[R, 20]) con elementos aleatorios entre 0 y 1. 
 b) Comprobar si son linealmente independientes 
 c) Determinar una base ortogonal del subespacio que generan, y comprobar el resultado.
 d) Determinar una base ortonormal del subespacio que generan, y comprobar el resultado.
 Sugerencias:  vec=Table[Table[Random[],{i,15}],{j,10}],
Det[vec.Transpose[vec]],  Y utilizar gs o Orthogonalize. 

5)  Encontrar la expresion general para la sucesión (Subscript[a, n + 1])=(Subscript[a, n])+2 (Subscript[a, n - 1

]) partiendo de (Subscript[a, 0]), (Subscript[a, 1]) genericos. "

a=.;v={a,b}

v[[2]]

a[[All,2]]

a={{1,1},{1,0}};Clear[n,a1,a0];{z,zv}=Eigensystem[a];s=Transpose[zv];si=Inverse[s];
an=Simplify[s.DiagonalMatrix[z^n].si.{a1,a0}][[2]]

(* caso usual de Fibonacci: *)

a0=0;a1=1;anf=Simplify[an]

ListPlot[Table[anf,{n,20}]]

tab=Table[Simplify[anf],{n,0,20,1}]

ListPlot[tab]

a={{1,2,3},{2,3,4},{5,-6,7}};

(* Sistemas Lineales  *)

Solve[{x+y+z==1,x-y-2z==2,x+y-2z==3},{x,y,z}]

Solve[{x+y+z==1,x-y-2z==2,x+y^2-2z==3},{x,y,z}]

Solve[{x+y+z==1,x-y-2z==2,x+5 y-2z==3},{x,y,z}]

(* Ecuaciones Diferenciales - resonancia (sin roce) *)

{yy[t_,w_]}=y[t]/.sol

t=.;y=.;Simplify[DSolve[y''[t]+9*y[t]==Cos[3 t],y[t],t]]

Clear[t,y,w];sol=Simplify[DSolve[{y''[t]+9*y[t]==Cos[w* t]},y[t],t]]

Clear[t,y,w];sol=Simplify[DSolve[{y''[t]+9*y[t]==Cos[w* t],y[0]==0,y'[0]==0},y[t],t]]

{yy[t_,w_]}=y[t]/.sol

ListAnimate[Table[Show[{Plot[yy[t,w],{t,0,30},PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],Graphics[{Text[StyleForm["w=",16],{3,2.5}],Text[StyleForm[w,16],{6,2.5}]}]}],{w,0.02,4.01,.25}]]

(* Resonancia con roce *)

Clear[t,y,w,ga];sol=Simplify[DSolve[{y''[t]+2*ga*y'[t]+4*y[t]==F*Cos[w* t]},y[t],t]]

Clear[t,y,w,ga];sol=Simplify[DSolve[{y''[t]+2*ga*y'[t]+4*y[t]==F*Cos[w* t],y[0]==0,y'[0]==0},y[t],t]]

{yy[t_,w_]}=y[t]/.sol

ga=.1;F=1;ListAnimate[Table[Show[{Plot[yy[t,w],{t,0,30},PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],Graphics[{Text[StyleForm["w=",16],{3,2.5}],Text[StyleForm[w,16],{6,2.5}]}]}],{w,0.02,4.01,.25}]]

w=2;F=1;ListAnimate[Table[Show[{Plot[ga=gaa;yy[t,w],{t,0,30},PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],Graphics[{Text[StyleForm["γ=",16],{3,2.5}],Text[StyleForm[ga,16],{6,2.5}]}]}],{gaa,0.02,4.01,.25}]]

yy[t,w]

ListAnimate[Table[Show[{Plot[yy[t,w],{t,0,30},PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],Graphics[{Text[StyleForm["w=",16],{3,2.5}],Text[StyleForm[w,16],{6,2.5}]}]}],{w,0.01,6,.25}]]

w=3.01;tabt=Table[{t=tt,yy[tt,w]},{tt,0,30,0.05}];ListAnimate[Table[Show[{ListPlot[Take[tabt,10*i],Joined->True,Filling->Axis,PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],
Graphics[{Text[StyleForm["w=",16],{3,2.5}],Text[StyleForm[w,16],{6,2.5}]}]}],{i,Length[tabt]/10}]]

Plot3D[yy[t,w],{t,0,18},{w,0,5},PlotRange->{{0,20},{0,5},{-3,3}},AxesLabel->{"t","w","y"}]

w=3.01;tabt=Table[{t=tt,yy[tt,w]},{tt,0,30,0.05}];ListAnimate[Table[Show[{ListPlot[Take[tabt,10*i],Joined->True,Filling->Axis,PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],
Graphics[{Text[StyleForm["w=",16],{3,2.5}],Text[StyleForm[w,16],{6,2.5}]}]}],{i,Length[tabt]/10}]]

w=3.01;

ListAnimate[Table[ListPlot[Table[{t=tt,yy[tt,w]},{tt,0,ttt,0.05}],Joined->True,PlotRange->{{0,30},{-3.5,3.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],{ttt,0.,30,.1}]]

Clear[t,y,w,ga];sol=Simplify[DSolve[{y''[t]+2*ga*y'[t]+9*y[t]==0},y[t],t]]

Clear[t,y,w,ga];sol=Simplify[DSolve[{y''[t]+2*ga*y'[t]+9*y[t]==0,y[0]==1,y'[0]==0},y[t],t]]

{yr[t_,ga_]}=y[t]/.sol;

yr[t,ga]

ListAnimate[Table[Show[{Plot[yr[t,ga],{t,0,10},PlotRange->{{0,10},{-1.5,1.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],Graphics[{Text[StyleForm["γ=",16],{3,1.2}],Text[StyleForm[ga,16],{3.85,1.2}]}]}],{ga,0.01,6.01,.5}]]

ga=.5;tabt=Table[{t=tt,yr[tt,ga]},{tt,0,10,0.025}];ListAnimate[Table[Show[{ListPlot[Take[tabt,10*i],Joined->True,Filling->Axis,PlotRange->{{0,10},{-1.5,1.5}},Frame->True,FrameLabel->{StyleForm["t",16],StyleForm["y(t)",16]},PlotStyle->Thick],
Graphics[{Text[StyleForm["γ=",16],{3,1.2}],Text[StyleForm[ga,16],{3.85,1.2}]}]}],{i,Length[tabt]/10}]]

(* Sistemas de ecuaciones diferenciales lineales *)

Clear[y,w,x,t];sol=Simplify[DSolve[{x'[t]==y[t],y'[t]==-x[t]+Cos[w*t]},{x[t],y[t]},t]]

(* con condiciones iniciales *)

Clear[y,w,x,t,f];sol=Simplify[DSolve[{x'[t]==y[t],y'[t]==-x[t]+f*Cos[w*t],y[0]==1,x[0]==0},{x[t],y[t]},t]]

{{xs[t_,w_],ys[t_,w_]}}={x[t],y[t]}/.sol

w=.8;f=1;Plot[{xs[t,w],ys[t,w]},{t,0,100},PlotStyle->Thick]

w=.8;ParametricPlot[{xs[t,w],ys[t,w]},{t,0,300},PlotStyle->Thick,PlotRange->{{-6,6},{-6,6}}]

Sound[Play[Sin[2*Pi*440*t],{t,0,6}]]

Sound[{SoundNote["C"],SoundNote["G"]}]

PENDULUM

y0=.;DSolve[{y''[t]==-w^2*Sin[y[t]],y[0]==0,y'[0]==v0},y[t],t]

LINEAR APPROXIMATION

y0=.;DSolve[{y''[t]==-w^2*y[t],y[0]==0,y'[0]==v0},y[t],t]

<<PhysicalConstants`

ProtonMass

ElectronMass