Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download
Project: WISB251
Views: 87
Kernel: Python 3 (Ubuntu Linux)

Inleveropdracht 3

In deze opdracht gaan we het randwaardeprobleem u(x)=f(x),voorx(0,2π) u''(x) = f(x), \quad \text{voor} \quad x \in (0,2\pi) en u(0)=u(2π)=1u(0) = u(2\pi) = 1 numeriek oplossen.

We discretiseren het interval hiertoe met N+1N+1 roosterpunten xn=nhx_n=n h, voor n{0,...,N}n\in\{ 0,...,N\} en met h=2π/Nh=2\pi /N en gaan op zoek naar een benadering van de oplossing, u~\widetilde{u}.

a) Laat zien dat de oplossing van het randwaardeprobleem voldoet het volgende stelsel vergelijkingen Au=f+e,A \mathbf{u}={\mathbf{f}} + \mathbf{e}, waar u{\mathbf{u}} een vector is met elementen u(xn)u(x_n) voor n=1,,N1n = 1, \ldots, N-1, f{\mathbf{f}} de bijbehorende waarden van ff bevat, e\mathbf{e} de truncatiefout is en A=N24π2(2112112112112).A= \frac{N^2}{4\pi^2} \left( \begin{matrix} -2 & 1 & & & & & \\ 1 & -2 & 1 & & & & \\ & 1 & -2 & 1 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & & 1 & -2 & 1 \\ & & & & & 1 & -2 \end{matrix} \right). Geef ook een uitdrukking voor de truncatiefout.

Hint: Beredeneer eerst welke benadering voor de tweede afgeleide u(x)u''(x) in ieder punt xnx_n gebruikt is. Denk ook aan het toepassen van de randvoorwaarden in x0x_0 en xNx_N.

antwoord: Met behulp van een Taylor expansie van u(xi±1)u(x_{i\pm 1}) rond xix_i vinden we u(xi)=(u(xi1)2u(xi)+u(xi1))/h2h2u(4)(ξi)/12.u''(x_i) = (u(x_{i-1}) -2u(x_i) + u(x_{i-1}))/h^2 - h^2 u^{(4)}(\xi_i)/12. In de eerste term herkennen we de coefficienten in te matrix. Gebruikmakend van het feit dat u(xi)=f(xi)u''(x_i) = f(x_i) vinden we dat e(xi)=h2u(4)(ξi)/12e(x_i) = h^2 u^{(4)}(\xi_i)/12.

1 punt

Het blijkt dat de eigenwaarden van AA zijn gegeven door λj=N2π2sin2(jπ2N)\lambda_j = -\frac{N^2}{\pi^2}\sin^2\left(\frac{j\pi}{2N}\right) voor j=1,,N1j = 1, \ldots , N-1. Je mag dit gebruiken bij het beantwoorden van de volgende vragen.

b) We kunnen nu de oplossing van het randvoorwaardeprobleem benaderen met u~\widetilde{\mathbf{u}} die voldoet aan Au~=fA\widetilde{\mathbf{u}}=\mathbf{f}. Geef een bovengrens voor de discretisatiefout uu~2\|\mathbf{u} - \widetilde{\mathbf{u}}\|_2 van de vorm uu~g(N) \|\mathbf{u} - \widetilde{\mathbf{u}}\| \leq g(N) en laat vervolgens zien dat de discretisatiefout van order N3/2N^{-3/2} is, dwz uu~=O(N3/2). \|\mathbf{u} - \widetilde{\mathbf{u}}\| = \mathcal{O}(N^{-3/2}).

Hint: Kijk naar het verschil A(uu~)A(\mathbf{u} - \widetilde{\mathbf{u}}) en gebruik wat je weet over de truncatiefout en de eigenwaarden van de matrix

antwoord: We kijken naar het verschil: A(uu~)=eA(\mathbf{u} - \widetilde{\mathbf{u}}) = \mathbf{e} en vinden dat uu~=A1e\mathbf{u} - \widetilde{\mathbf{u}} = A^{-1}\mathbf{e}. Gebruikmakend van de eigenschappen van de norm krijgen we de volgende afschatting uu~A1e.\|\mathbf{u} - \widetilde{\mathbf{u}}\| \leq \|A^{-1}\|\|\mathbf{e}\|. (1 punt)

Eerst gaan we een uitdrukking voor A1\|A^{-1}\| vinden. Als we de 2-norm gebruiken en de relatie B2=ρ(BTB)\|B\|_2 = \sqrt{\rho(B^TB)} (waarbij ρ\rho de spectrale radius is) krijgen we A12=λmax(A1)=1/λmin(A)\|A^{-1}\|_2=|\lambda_{\max}(A^{-1})| = |1/\lambda_{\min}(A)| (omdat AA symmetrisch is). De eigenwaarden van AA zijn gegeven en we vinden A12=π2N2sin2(jπ2N)\|A^{-1}\|_2 = \frac{\pi^2}{N^2\sin^2\left(\frac{j\pi}{2N}\right)} (1/2 punt)

Om de tweede term af te schatten maken we gebruik van de ongelijkheid e2N1e\|\mathbf{e}\|_2 \leq \sqrt{N-1}\|\mathbf{e}\|_{\infty}. Hieruit krijgen we e2(2π)2N112N2M,\|\mathbf{e}\|_2 \leq (2\pi)^2\frac{\sqrt{N-1}}{12 N^2}M, waarbij M=maxiu(4)(xi).M = \max_{i} |u^{(4)}(x_i)|. (Als je gebruikt dat e2e1(N1)e\|\mathbf{e}\|_2 \leq \|\mathbf{e}\|_1 \leq (N-1)\|\mathbf{e}\|_{\infty} en op een factor (N1)(N-1) uitkomt is dat ook goed.) (1/2 punt)

We gaan nu de orde van grootte van de fout bepalen. We gebruiken hier dat limsin(x)x=1\lim_{\rightarrow\infty} \frac{\sin(x)}{x} = 1 en vinden zo het gewenste antwoord. (Als je op O(N1)\mathcal{O}(N^{-1}) uitkomt is dat ook goed).

(1/2 punt)

c) We gaan nu de oplossing van het stelsel vergelijkingen benaderen door het stelset iteratief op te lossen met een vastepuntiteratie u~(m+1)=u~(m)+α(fAu~(m)). \widetilde{\mathbf{u}}^{(m+1)} = \widetilde{\mathbf{u}}^{(m)} + \alpha \left({\mathbf{f}} - A\widetilde{\mathbf{u}}^{(m)}\right). Laat zien dat de iteratieconvergeert voor α(2π2N2sin2(π(N1)/(2N)),0)\alpha \in (-\frac{2\pi^2}{N^2\sin^2(\pi(N-1)/(2N))},0) en geef ook een optimale waarde voor α\alpha. Hint: kijk nog eens naar opgave 7.7.9

antwoord: We krijgen de volgende uitdrukking voor de fout u~(m+1)u~=(IαA)(u~(m+1)u~)=(IαA)m+1(u~(0)u~).\widetilde{\mathbf{u}}^{(m+1)} - \widetilde{\mathbf{u}} = \left(I - \alpha A\right)\left(\widetilde{\mathbf{u}}^{(m+1)} - \widetilde{\mathbf{u}}\right) = \left(I - \alpha A\right)^{m+1}\left(\widetilde{\mathbf{u}}^{(0)} - \widetilde{\mathbf{u}}\right). Deze iteratie convergeert als ρ(IαA)<1\rho(I - \alpha A) < 1 (zie pagina 180 van het boek). De eigenwaarden van IαAI-\alpha A zijn 1αλj1 - \alpha\lambda_j. De grootste (in absolute zin) is 1αλN11 - \alpha\lambda_{N-1}. De grenzen zijn dus α=0\alpha = 0 en α=2/λN1\alpha = 2/\lambda_{N-1}. Invullen geeft het gewenste antwoord. (1/2 punt).

De optimale α\alpha geeft een zo klein mogelijke spectrale radius. Omdat IαAI-\alpha A zowel positieve als negatieve eigenwaarden kan hebben verkrijgen we de kleinst mogelijke als de meest negatieve gelijk is is magnitude als de meest positieve, ofwel 1αλ1=(1αλN1).1 - \alpha\lambda_1 = - (1-\alpha\lambda_{N_1}). Hieruit volgt α=2λ1+λN1.\alpha = \frac{2}{\lambda_1 + \lambda_{N-1}}. (1/2 punt).

d) Laat zien dat de convergentiefout na MM iteraties is begrensd door u~(M)u~2ρMu~2, \|\widetilde{\mathbf{u}}^{(M)} - \widetilde{\mathbf{u}}\|_2 \leq \rho^M\|\widetilde{\mathbf{u}}\|_2, ervan uitgaande dat u~(0)=0\widetilde{\mathbf{u}}^{(0)} = \mathbf{0} en geef een uitdrukking voor ρ\rho.

antwoord: Uit het antwoord in c) krijgen we direct

u~(M)u~2IαA2Mu~2,\|\widetilde{\mathbf{u}}^{(M)} - \widetilde{\mathbf{u}}\|_2 \leq \|I-\alpha A\|_2^{M}\|\widetilde{\mathbf{u}}\|_2,

De spectrale radius van IαAI-\alpha A is (met de optimale α\alpha) ρ=12λN1λ1+λN1=λ1λN1λ1+λN1.\rho = 1 - 2\frac{\lambda_{N-1}}{\lambda_1 + \lambda_{N-1}} = \frac{\lambda_1 - \lambda_{N-1}}{\lambda_1 + \lambda_{N-1}}. (1/2 punt)

e) Geef nu een bovengrens voor de totale fout uu~(M)\|\mathbf{u} - \widetilde{\mathbf{u}}^{(M)}\| als functie van MM and NN en geef aan hoe de fout zich asymptotisch gedraagt als M,NM,N \rightarrow \infty.

antwoord: We herschrijven de fout eerst als (uu~)(u~(M)u~)2\|(\mathbf{u} - \widetilde{\mathbf{u}}) - (\widetilde{\mathbf{u}}^{(M)} - \widetilde{\mathbf{u}})\|_2 en gebruiken de driehoeksongelijkheid: uu~(M)2uu~2+u~(M)u~2.\|\mathbf{u} - \widetilde{\mathbf{u}}^{(M)}\|_2 \leq \|\mathbf{u} - \widetilde{\mathbf{u}}\|_2 + \|\widetilde{\mathbf{u}}^{(M)} - \widetilde{\mathbf{u}}\|_2. (1/2 punt)

Van de eerste term weten we dat deze O(N3/2)\mathcal{O}(N^{3/2}) is. Van de tweede term merken we het volgende op: Als NN\rightarrow\infty voor vaste MM, dan ρ=O(1)\rho = \mathcal{O}(1). Voor vaste NN hebben we dat ρ<1\rho < 1, dus ρM=O(1)\rho^{M} = \mathcal{O}(1) als MM\rightarrow\infty. Als laatste merken we op dat u~2A12f2N1A12f\|\widetilde{\mathbf{u}}\|_2 \leq \|A^{-1}\|_2\|\mathbf{f}\|_2 \leq \sqrt{N-1}\|A^{-1}\|_2\|\mathbf{f}\|_\infty wat we kunnen afschatten op O(N1/2)\mathcal{O}(N^{1/2}). (1/2 punt)

Om er voor te zorgen dat de tweede term kleiner wordt als NN\rightarrow\infty zal MM dus ook moeten groeien. (Ik verwacht niet dat iemand zo ver is gekomen; ik heb zelf ook nog niet uitgewerkt hoe snel MM moet groeien om die hele fout te laten krimpen als NN\rightarrow\infty....)

f) Ga nu door middel van een numeriek experiment na hoe goed de bovengrens is die je bij e) hebt afgeleid. Neem hiervoor f(x)=sin(x)f(x) = -\sin(x), zodat de echte oplossing is gegeven door u(x)=sin(x)u(x) = \sin(x). Laat een grafiek zien van de discretisatiefout, convergentiefout en totale fout als een functie van NN voor verschillende MM. Bespreek wat je observeert in de grafieken.

Hieronder vind de benodigde code om je op weg te helpen.

import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as splinalg import matplotlib.pyplot as plt N = 20 M = 200 # rooster h = 2*np.pi/N xi = h*np.linspace(1,N-1,N-1).reshape((N-1,1)) # matrix A = (1/h**2)*sp.diags([1, -2, 1], [-1, 0, 1], shape=(N-1,N-1)) # bronterm f = -np.sin(xi) # echte oplossing u = np.sin(xi) # echte oplossing van het stelsel ut = splinalg.spsolve(A,f).reshape((N-1,1)) # iteratieve oplossing alpha = -2*np.pi**2/N**2 ui = np.zeros((N-1,1)) for k in range(M): ui = ui + alpha*(f - A@ui) # fout print('discretisatiefout : ', np.linalg.norm(ut - u)) print('convergentiefout : ', np.linalg.norm(ui - ut)) print('total fout : ', np.linalg.norm(ui - u)) # grafiek van oplossingen plt.plot(xi,u,xi,ut,xi,ui) plt.legend(('sin(x)','ut','ui'))
discretisatiefout : 0.0261375434243 convergentiefout : 0.000139585395153 total fout : 0.0259979580291
/usr/local/lib/python3.5/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format SparseEfficiencyWarning)
<matplotlib.legend.Legend at 0x7fb577809f98>
Image in a Jupyter notebook
def error(N,M=1): # rooster h = 2*np.pi/N xi = h*np.linspace(1,N-1,N-1).reshape((N-1,1)) # matrix A = (1/h**2)*sp.diags([1, -2, 1], [-1, 0, 1], shape=(N-1,N-1)) # bronterm f = -np.sin(xi) # echte oplossing u = np.sin(xi) # echte oplossing van het stelsel ut = splinalg.spsolve(A,f).reshape((N-1,1)) # iteratieve oplossing alpha = -2*np.pi**2/N**2 ui = np.zeros((N-1,1)) for k in range(M): ui = ui + alpha*(f - A@ui) return (np.linalg.norm(u - ut), np.linalg.norm(ui - ut), np.linalg.norm(ui - u))

antwoord: Hieronder een plot van de discretisatiefout als functie van NN. We observeren inderdaad dat de fout afneemt als N3/2N^{-3/2}. (1 punt)

# discretization error ns = np.linspace(10,100,91) e = [] for n in ns: en = error(int(n),1) e.append(en[0]) plt.loglog(ns,e,ns,ns**(-3/2)) plt.legend(('d. error','bound'))
/usr/local/lib/python3.5/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format SparseEfficiencyWarning)
<matplotlib.legend.Legend at 0x7fb577812588>
Image in a Jupyter notebook

antwoord: Hieronder de convergentiefout waarin we lineaire convergentie observeren, zoals verwacht. De bovengrens is wel erg pessismistisch, in praktijk neemt de fout sneller af. (1 punt)

# approximation error ns = 10 ms = np.linspace(0,100,100) e = [] for m in ms: en = error(int(ns),int(m)) e.append(en[1]) plt.semilogy(ms,e,ms,((1 - (np.pi/2/ns)**2)**ms)*e[0]) plt.legend(('a error','bound'))
/usr/local/lib/python3.5/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format SparseEfficiencyWarning)
<matplotlib.legend.Legend at 0x7fb5776e22e8>
Image in a Jupyter notebook

antwoord: Hieronder een grafiek van de convergentiefout als functie van NN. We zien dat de fout langzamer afneemt naarmate NN toeneemt. Bijzonder is ook dat de fout toe begint te nemen naarmate MM toeneemt. Duidelijk te zien is het O(1)\mathcal{O}(1) gedrag van de totale fout als MM\rightarrow\infty. (1 point)

# total error ms = np.linspace(0,1000,100) ns = [10,20,40] for n in ns: e = [] for m in ms: en = error(int(n),int(m)) e.append(en[2]) plt.semilogy(ms,e) plt.legend(ns)
/usr/local/lib/python3.5/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format SparseEfficiencyWarning)
<matplotlib.legend.Legend at 0x7fb577754a20>
Image in a Jupyter notebook

Tot slot vragen we ons af hoe snel MM moet groeien met NN als we een convergende totale fout willen. In de grafiek hieronder zien we dat M=O(N2)M = \mathcal{O}(N^2) volstaat om de totale fout met O(N3/2)\mathcal{O}(N^{-3/2}) te laten afnemen. (1 punt)

# total error ns = np.linspace(10,100,91) e = [] for n in ns: en = error(int(n),int(n**(2))) e.append(en[2]) plt.loglog(ns,e,ns,ns**(-3/2)) plt.legend(('total error','O(N^(-3/2))'))
/usr/local/lib/python3.5/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format SparseEfficiencyWarning)
<matplotlib.legend.Legend at 0x7fb5772e9c18>
Image in a Jupyter notebook