(* ::Package:: *)

(*
Rechnernutzung in der Physik, WS23/24

Nicht-linearer Oszillator.
*)

(* ***
Get["kap_3_non-lin-osz.m"];
NichtLinOszil[4];
?u
?b

TrigReduce[Cos[t]^2]
TrigReduce[Cos[t]*Cos[2*t]]

ex1 = Series[Cos[x],{x,0,3}]
ex2 = Series[Sin[x],{x,0,5}]
ex1 * ex2
*** *)


(* n: verwende x^n-Term in V(x) *)
NichtLinOszil[n_] := Module[
    {V,R,Eq,i,Et},
(* Potential: *)
    Clear[x,w2,b,u];
    V = Series[Sum[c[i]/i*x^i,{i,3,n}],{x,0,n}];
    Print["V = ",V];
(* rechte Seite der Differentialgleichung: *)
    R = -D[V,x];
    Print["R = ",R];
(* Reihenentwicklung der "Frequenz"; nur gerade Terme in a: *)
    w2 = Series[1+Sum[u[i]*a^i,{i,2,n,2}],{a,0,n-2}];
(* Ansatz fuer die allgemeine Loesung: Zur Ordnung a^i tauchen Terme Cos[j*t] fuer j<=i auf.
   j ist gerade/ungerade falls i gerade/ungerade ist.
   Der Fall j=1 taucht nicht auf, da er per Definition im fuehrenden Term enthalten ist. *)
    x[t] = Series[ a*Cos[t] + Sum[ a^i*Sum[b[i,j]*Cos[j*t] ,{j, If[EvenQ[i],0,3],i,2} ], {i,2,n-1}], {a,0,n-1}];
(* Bewegungsgleichung: *)
    Eq = Map[TrigReduce, ExpandAll[w2*D[x[t],{t,2}]+x[t]-(R/.x->x[t])] ];
(* Alle Koeffizienten von a^i, i=2,...,n-1 muessen Null sein.
   Falls i ungerade ist, gibt es Cos[t]-Terme. Die Bedingung, dass der Vorfaktor verschwindet
   liefert Korrekturen zur Frequenz. Alle anderen Terme liefern Koeffizienten fuer x(t). *)
    Do[   Eqi = SeriesCoefficient[Eq,i] /. Cos[j_.*t]->z^j;
          (*Print["i=",i,": Eqi=",Eqi];*)
  		If[OddQ[i], u[i-1] = u[i-1] /. Solve[ Coefficient[Eqi,z,1]==0, u[i-1]][[1]] ];
	      Do[ b[i,j] = b[i,j] /. Solve[Coefficient[Eqi,z,j]==0,b[i,j]][[1]], {j,If[EvenQ[i],0,3],i,2} ]
	,{i,2,n-1} ];
(* "Frequenz": *)
    w2 = w2;
    Print["omega^2 = ",w2];
(* x(t): *)
    xres = Collect[x[t],{a,Cos[xx_]},Expand];
    Print["x[t] = ",xres];
(* Check: Energieerhaltung: *)
    Et = Map[TrigReduce,ExpandAll[ w2*D[x[t],t]^2/2 + x[t]^2/2 + (V/.x->x[t])]];
    Print["Energie = ",Et];

];