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

Poissonvergelijking

De poissonvergelijking op [0,1]2[0,1]^2 is gegeven door

2u(x,y)=f(x,y),\nabla^2 u(x,y) = f(x,y),

waarbij ff een gegeven functie is. Om de oplossing uniek te maken leggen we ook randvoorwaarden op, bijvoorbeeld u(0,x)=u(1,x)=u(x,0)=u(x,1)=0u(0,x) = u(1,x) = u(x,0) = u(x,1) = 0. Deze vergelijking komt veel voor, en kan bijvoorbeeld worden gebruikt om de warmteverdeling in een kamer als gevolg van een hittebron ff te beschrijven. Met f=0f = 0 en andere randvoorwaarden komt dit probleem in de aerodynamica voor.

We gaan deze vergelijking numeriek oplossen door de Laplace operator te benaderen:

2u(x,y)(u(x+h,y)+u(x,y+h)4u(x,y)+u(xh,y)+u(x,yh))/h2,\nabla^2 u(x,y) \approx \left(u(x+h,y) + u(x,y+h) - 4u(x,y) + u(x-h,y) + u(x,y-h)\right)/h^2,

waar hh de stapgrootte is. Als we nu uiju(ih,jh)u_{ij} \equiv u(ih, jh) invoeren kunnen we Poissonvergelijking als een stelsel vergelijkingen schrijven

Lu=f.L\mathbf{u} = \mathbf{f}.

Omdat u=0u = 0 op de rand kunnen we de termen op de rand negeren. We krijgen nu een stelsel van n2n^2 vergelijkingen met n2n^2 onbekenden, met h=1/(n+1)h = 1/(n+1).

Omdat deze matrix ijl is (er zijn slechts 5 niet-nullen in elke rij) is het niet erg aantrekkelijk om dit stelsel op te lossen met een eliminatie procedure (waarom niet?). We gaan daarom kijken naar iteratieve methoden. De eenvoudigste is de volgende vastepunt iteratie

un+1=un+ω(fLun),\mathbf{u}_{n+1} = \mathbf{u}_{n} + \omega \left(\mathbf{f} - L\mathbf{u}_n\right),

waar ω\omega de relaxatieparameter is. We hebben eerder gezien dat de convergentiesnelheid afhangt van de spectrale radius van de matrix Iω1LI - \omega^{-1}L. We zullen echter zien dat dit niet alleszeggend is; de convergentiesnelheid blijkt ook af te hangen van de gegeven f\mathbf{f}, dat wil zeggen, voor sommige f\mathbf{f} is de convergentie veel sneller dan de spectrale radius doet vermoeden.

Het probleem zit 'm in dit geval in de grote spreiding van de eigenwaarden (i.e., het verschil tussen de grootste en de kleinste eigenwaarden). Om hier iets tegen te doen gaan we kijken naar iterative methoden van de vorm

un+1=un+M1(fLun),\mathbf{u}_{n+1} = \mathbf{u}_{n} + M^{-1}\left(\mathbf{f} - L\mathbf{u}_n\right),

waar we de matrix MM zo kiezen dat M1LIM^{-1}L \approx I. Een aantal mogelijk keuzes worden in hoofdstuk 7.2 van het boek besproken.

Vragen voor het verslag

Basis

  • Geef een uitdrukking voor de discretisatiefout

  • Beschrijf hoe de matrix LL eruit ziet.

  • Beschrijf wat er gebeurt als je het stelsel zou oplossen met een directe (eliminatie) methode.

  • De eigenwaarden van LL zijn gegeven door λij=4h2(sin2(iπh/2)+sin2(jπh/2))\lambda_{ij} = -\frac{4}{h^2}\left(\sin^2(i\pi h/2) + \sin^2(j\pi h/2)\right) voor i,j=1,2,,ni,j = 1, 2, \ldots,n, bereken de eigenwaarden ook numeriek om dit te controleren. Hoe zien de bijbehorende eigenvectoren eruit? Gebruik numpy.linalg.eigvals om deze uit te rekenen.

  • Wat is de optimale ω\omega voor een zo snel mogelijke convergentie?

  • Hoe hangt de convergentie af van hh in het ergste geval?

Extra

  • Probeer de standaard iterative methode (met de optimale ω\omega) eens op f=1\mathbf{f} = \mathbf{1} en een f\mathbf{f} met willekeurige elementen (f = randn(n^2,1)). Wat kun je zeggen over de convergentie, hoe is dit vergeleken met de theoretische voorspelling?

  • Leg uit waarom dit gebeurt. Hint: kijk eens naar de eigenwaarden en bijbehorde eigenvectoren van LL.

  • Probeer verschillende keuzes voor MM en vergelijk de convergentie. Bespreek voor- en nadelen van ieder.

import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla import matplotlib.pyplot as plt n = 100 h = 1/(n+1) # Laplace operator D = (1/h**2)*sp.diags([1, -2, 1],[-1, 0, 1],(n,n)) I = sp.diags([1],[0],(n,n)) L = sp.kron(D,I) + sp.kron(I,D) # bekijk hoe de matrices eruit zien plt.spy(L)
<matplotlib.lines.Line2D at 0x7f3d01786a58>
Image in a Jupyter notebook
L
<9x9 sparse matrix of type '<class 'numpy.float64'>' with 33 stored elements in Compressed Sparse Row format>
# f = np.random.randn(n**2).reshape((n**2,1)) # vastepuntiteratie N = 100 omega = -1e-4 u = np.zeros((n**2,1)) r = np.zeros(N) for k in range(N): r[k] = np.linalg.norm(L@u - f) u = u + omega*(f - L@u) # plot fig1, axs1 = plt.subplots(1,2) axs1[0].loglog(r) axs1[1].imshow(u.reshape((n,n)))
<matplotlib.image.AxesImage at 0x7f3d01f0a0b8>
Image in a Jupyter notebook
# grid xx,yy = np.meshgrid(np.linspace(h,1-h,n),np.linspace(h,1-h,n)) # bronterm f = -(2*np.pi**2)*(np.sin(np.pi*xx)*np.sin(np.pi*yy)) # echte oplossing u = (np.sin(np.pi*xx)*np.sin(np.pi*yy)) # numerieke benadering ut = spla.spsolve(L,f.reshape(n**2,1)) # plot fig2, axs2 = plt.subplots(1,3) axs2[0].imshow(u,vmin=0,vmax=1) axs2[0].set_title('echte oplossing') axs2[1].imshow(ut.reshape((n,n)),vmin=0,vmax=1) axs2[1].set_title('numerieke oplossing') axs2[2].imshow(ut.reshape((n,n)) - u,vmin=-1,vmax=1) axs2[2].set_title('fout')
Text(0.5,1,'fout')
Image in a Jupyter notebook