Exit[]

CurrentValue[$FrontEnd,CellContext]=Notebook;

Zkuste porozumět násldeujcímu postupu řešení dif. rovnice

Nejprve rovnice, kterou řešíme

rovnice = D[f[x],x,x]==-Cos[Pi x] f[x]
okrajPodminky = f[-1]==1 && f[1]==1

Kupodivu, rovnice má přesné řešení ve třídě velmi speciálních funkcí MathieuC
To se hodí při zkoumání přesnosti numerického řešení

DSolve[ rovnice && okrajPodminky,f[x],x]
reseni =f[x]/.%[[1]]
Plot[ reseni,{x,-1,1},PlotTheme->"Detailed"]

Nyní budeme hledat aproximaci řešení s použitím polynomů
Teorie říká, že máme vzít Čebyševovy polynomy 
Jak moc je konkrétní polynom zastoupen v řešení představuje neznámé (Subscript[a, 0])... (Subscript[a, n])

Problém 2.1. Namalujte do jednoho grafu na intervalu (-1,1) polynomy (Subscript[T, 15])(x) a (x^15)

nMax =26;
obecnyPolynom = Sum[a[k]ChebyshevT[k,x],{k,0,nMax}]
nezname = Table[a[k],{k,0,nMax}]

Dosazením polynomu do rovnice získáme rovnici, kterou by měl splňovat
Polynom musí splňovat i okrajové podmínky

residual = rovnice /. f->Function[x,Evaluate@obecnyPolynom]
okrPodminkyProa = okrajPodminky /. f->Function[x,Evaluate@obecnyPolynom]

Máme 2 okrajové podmínky, ale n+1 neznámých (Subscript[a, 0])... (Subscript[a, n])
Potřebujeme n+1 rovnic, tedy n-1 vyrobíme dosazením vhodných x do residualu
Těmto bodům se říká kolokační a celé metodě se říká pseudospektrální metoda řešení ODE

Jakmile máme  dost rovnic, vyřešíme je

pocetNeznamych == Length[nezname]
nC = nMax-1;

colocationPoints = Table[ Cos[  (2k+1) Pi/2 / nC ],{k,0,nC-1}]//N
rovnicProOkrajovePodminky == Length @ okrPodminkyProa
pocetBoduX == Length[colocationPoints]

(* dame dohromady vsechny rovnice *)
residual/. ({x->#}&)/@ colocationPoints;
(%/.List -> And) && okrPodminkyProa;
(* a vyresime je *)
numReseniA = Solve[%,nezname]

Protože je přesné řešení sudá funkce, vycházejí lichá a hodně malá, lze se jich zbavit funkcí Chop

fNum = obecnyPolynom/.numReseniA//Chop

Můžeme nakreslit přesné a přibližné řešení a také nalézt odhad maximální odchylky

Plot[ {fNum,reseni},{x,-1,1},PlotTheme->"Detailed",PlotLegends->{"približně","přesně"}]
MaximalniChyba[n==nMax] == (Table[fNum-reseni,{x,-1,1,1/1000}]//Abs//Max)

Pseudospektrální metoda za příznivých okolností vykazuje chybu srovnatelnou s velikostí posledních  členů rozvoje. Viz též graf Čebyševových polynomů.

posledniClenyRozvoje == {a[nMax-2],a[nMax-1],a[nMax]}/.numReseniA

Problém 2.2. Namalujte s pomocí funkce ListLogPlot závislost | (Subscript[a, k])|  na indexu k pro dosatečně velké nMax tak jako je na grafu níže. Takovýto graf říká kolik členů rozvoje zhruba potřebujeme, abychom rovnici vyřešili s danou přesností.

ListLogPlot [ Abs@nezname/.numReseniA[[1]],Frame->True,GridLines->Automatic,PlotRange->Full,FrameLabel->{k,"|(Subscript[a, k])|"},BaseStyle->14]

Abs@nezname/.numReseniA[[1]]