Exit[]

CurrentValue[$FrontEnd,CellContext]=Notebook;

Demonstrace (pseudo)spektrální metody řešení ODR

( Vynikající kniha: https://depts.washington.edu/ph506/Boyd.pdf )

Funguje samozřejmě i pro regulární rovnice, zde je ovšem záměrně použita singulární ODR jako příklad, kde NDSolve selže. 
(Navíc pro singulární rovnice, co máme je druhá okrajová podmínka nadbytečná až škodlivá)

Nejprve (Besselova) diferenciální rovnice, kterou řešíme.

rovnice = x^2 f''[x]+x f'[x]+9x^2 f[x]==0
okrajPodminky = f[-1]==1 && f[1]==1

Rovnice má přesné řešení. To se hodí při zkoumání přesnosti numerického řešení

Problém 1.1  Namalujte toto přesné řešení na intervalu (-1,1) + schovejte výsledek do proměnné reseni.

DSolve[⋯,⋯]
reseni =⋯
Plot[ reseni,{x,-1,1},PlotTheme->"Detailed"]

Rovnice má je nicméně singulární a to brání použití NDSolve s BVP

NDSolve[ rovnice && okrajPodminky,f[x],{x,-1,1}]

To proto, že i pokud problém přeformulujeme jako  IVP, nedostaneme se přes x==0

NDSolve[ rovnice && f[-1]==1 && f'[-1]==-3.9114,f[x],{x,-1,1}]
GraphicsRow[ {Plot[ f[x]/.%,{x,-1,0}], Plot[ f[x]/.%,{x,-0.001,0}] }]

Jako numerickou alternativu použijeme pseudospektrální metodu.

Tedy budeme hledat aproximaci řešení v podobě polynomů
Teorie říká, že máme vzít Čebyševovy polynomy (Subscript[T, n])(x) nikoli (x^k)
Jak moc je konkrétní polynom zastoupen v řešení představuje neznámé koeficienty (Subscript[a, 0])... (Subscript[a, n]) 

Pozn. Správně bychom měli zabránit tomu, aby byly (Subscript[T, n]) převedeny do podoby polynomů v x

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

Problém 2.2. Demonstrujte potíže se zaokrouhlováním (konečnou přesností) na ChebyshevT[50,0.99]

Definujeme polynom s obecnými koeficinety. Ty budeme muset určit, proto je budeme chápat jako neznámé.

nMax =7;
presnost = 45;
obecnyPolynom = Sum[a[k]T[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 /. ⋯
okrPodminkyProa = okrajPodminky /.⋯

T a U tvoří uzavřený systém vzhledem k derivaci (jsou to siny a kosiny takže to nepřekvapí)

{ T'[x]==D[ ChebyshevT[k,x],x], T''[x]==D[ ChebyshevT[k,x],x,x]}

Problém 2.3. Zařiďte, aby se pravidla výše aplikovala na residual (definujet pleaseNoT)

bezderivaci:={⋯}

residual/.bezderivaci

Problém 2.4. Definujte numerické hodnoty T a U (pro celé k a číselné x) viz _?NumberQ

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í ODR

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

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

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

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

Problém 2.5. Namalujte polohu kolokačních bodů na ose x pro n=7, 8, 30

Definujeme funkci reprezentujici nalezene reseni

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

fNum =  obecnyPolynom/.numReseniA

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)

Problém 2.6. Namalujte rozdíl mezi přesným a přibližným řešením

Problém 2.5. Namalujte s pomocí funkce ListLogPlot závislost | (Subscript[a, k])|  na indexu k pro dosatečně velké nMax≈31 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í.