Minimalne ploskve (Laplaceova enačba)

Minimalne ploskve (Laplaceova enačba)

Naloga

Žično zanko s pravokotnim tlorisom potopimo v milnico, tako da se nanjo napne milna opna.

Radi bi poiskali obliko milne opne, razpete na žični zanki. Malo brskanja po fizikalnih knjigah in internetu hitro razkrije, da ploskve, ki tako nastanejo, sodijo med minimalne ploskve, ki so burile domišljijo mnogim matematikom in nematematikom. Minimalne ploskve so navdihovale tudi umetnike npr. znanega arhitekta Otto Frei, ki je sodeloval pri zasnovi Muenchenskega olimpijskega stadiona, kjer ima streha obliko minimalne ploskve.

Muenchenski Olimpijski park (slike so vzete iz wikipedie) Slika wikipedia

Matematično ozadje

Ploskev lahko predstavimo s funkcijo dveh spremenljivk $u(x,y)$, ki predstavlja višino ploskve nad točko $(x,y)$. Naša naloga bo poiskati funkcijo $u(x,y)$ na tlorisu žične mreže.

Funkcija $u(x,y)$, ki opisuje milno opno, zadošča matematična enačbi, znani pod imenom Poissonova enačba

\[\Delta u(x,y)=\rho(x,y).\label{eq:Poisson}\]

Funkcija $\rho(x,y)$ je sorazmerna tlačni razliki med zunanjo in notranjo površino milne opne. Tlačna razlika je lahko posledica višjega tlaka v notranjosti milnega mehurčka ali pa teže, v primeru opne, napete na žični zanki. V primeru minimalnih ploskev pa tlačno razliko kar zanemarimo in dobimo Laplaceovo enačbo

\[\Delta u(x,y) = 0\label{Laplace}.\]

Če predpostavimo, da je oblika na robu območja določena z obliko zanke, rešujemo robni problem za Laplaceovo enačbo. Predpostavimo, da je območje pravokotnik $[a, b]\times[c, d]$. Poleg Laplacove enačbe, veljajo za vrednosti funkcije $u(x, y)$ tudi robni pogoji:

\[\begin{array}{ll} u(x, c) = f_s(x) & u(x, d) = f_z(x)\\ u(a, y) = f_l(y) & u(b, y) = f_d(y),\\ \end{array}\]

kjer so $f_s, f_z, f_l$ in $f_d$ dane funkcije. Rešitev robnega problema je tako odvisna od območja, kot tudi od robnih pogojev.

Za numerično rešitev Laplaceove enačbe za minimalno ploskev dobimo navdih pri arhitektu Frei Otto, ki je minimalne ploskve raziskoval tudi z elastičnimi tkaninami.

Diskretizacija in linearni sistem enačb

Problema se bomo lotili numerično, zato bomo vrednosti $u(x,y)$ poiskali le v končno mnogo točkah: problem bomo diskretizirali. Za diskretizacijo je najpreprosteje uporabiti enakomerno razporejeno pravokotno mrežo točk. Točke na mreži imenujemo vozlišča. Zaradi enostavnosti bomo obravnavali le mreže z enakim razmikom v obeh koordinatnih smereh. Omejimo se le na pravokotna območja v ravnini $[a, b]\times[c, d]$. Interval $[a, b]$ razdelimo na $n+1$ delov, interval $[c, d]$ pa na $m+1$ delov in dobimo zaporedje koordinat

\[\begin{array}{lcl} a=x_0, x_1,&\ldots&x_{n+1}=b\\ c=y_0, y_1,&\ldots&y_{m+1}=d,\\ \end{array}\]

ki definirajo pravokotno mrežo točk $(x_i, y_j)$. Namesto funkcije $u: [a,b]\times[c,d]\to \mathbb{R}$ tako iščemo le vrednosti

\[u_{ij} = u(x_i, y_j),\quad i=1,\ldots n, j=1,\ldots m\]

sosednja vozlišča

Iščemo torej enačbe, ki jim zadoščajo elementi matrike $u_{ij}$. Laplaceovo enačbo lahko diskretiziramo z končnimi diferencami, lahko pa izpeljemo enačbe, če si ploskev predstavljamo kot elastično tkanino, ki je fina kvadratna mreža iz elastičnih nitk. Vsako vozlišče v mreži je povezano s 4 sosednjimi vozlišči. Vozlišče bo v ravnovesju, ko bo vsota vseh sil nanj enaka 0. Predpostavimo, da so vozlišča povezana z idealnimi vzmetmi in je sila sorazmerna z razliko. Če zapišemo enačbo za komponente sile v smeri $z$, dobimo za točko $(x_i, y_j, u_{ij})$ enačbo

\[u_{i-1,j}+u_{i,j-1} - 4u_{ij} + u_{i+1,j}+u_{i,j+1} = 0.\]

Za $u_{ij}$ imamo tako sistem linearnih enačb. Ker pa so vrednotsi na robu določene z robnimi pogoji, moramo elemente $u_{0j}$, $u_{n+1,j}$, $u_{i0}$ in $u_{im+1}$ prestaviti na desno stran in jih upoštevati kot konstante.

Matrika sistema linearnih enačb

Sisteme linearnih enačb običajno zapišemo v matrični obliki

\[A\mathbf{x}=\mathbf{b},\]

kjer je $A$ kvadratna matrika, $\mathbf{x}$ in $\mathbf{b}$ pa vektorja. Spremenljivke $u_{ij}$ razvrstimo po stolpcih v vektor.

Razvrstitev po stolpih

Eden od načinov, kako lahko elemente matrike razvrstimo v vektor, je, da stolpce matrike enega za drugim postavimo v vektor. Indeks v vektorju $k$ lahko izrazimo z indeksi $i,j$ v matriki s formulo $k = i+(n-1)j.$

Za $n=m=3$ dobimo $9\times 9$ matriko

\[L = \begin{bmatrix} -4& 1& 0& 1& 0& 0& 0& 0& 0\\ 1& -4& 1& 0& 1& 0& 0& 0& 0\\ 0& 1& -4& 0& 0& 1& 0& 0& 0\\ 1& 0& 0& -4& 1& 0& 1& 0& 0\\ 0& 1& 0& 1& -4& 1& 0& 1& 0\\ 0& 0& 1& 0& 1& -4& 0& 0& 1\\ 0& 0& 0& 1& 0& 0& -4& 1& 0\\ 0& 0& 0& 0& 1& 0& 1& -4& 1\\ 0& 0& 0& 0& 0& 1& 0& 1& -4\\ \end{bmatrix},\]

ki je sestavljena iz $3\times 3$ blokov

\[\begin{bmatrix}-4&1&0\cr 1&-4&1\cr 0&1&-4\end{bmatrix},\quad \begin{bmatrix}1&0&0\cr 0&1&0\cr 0&0&1\end{bmatrix}.\]

desne strani pa so

\[\mathbf{b} = -[u_{01}+u_{10}, u_{20}, \ldots u_{n0}+u_{n+1,1},u_{02}, 0,\ldots u_{n+1,2}, u_{03}, 0\ldots u_{n, m+1},u_{n,m+1}+u_{n+1,m}]^T.\]

Primer

Z, x, y = resi_robni_problem(x->sin(x), y->-sin(y),
            x->sin(x), y->-sin(y), 0.1, ((0, pi), (0, pi)))
surface(x, y, Z)
savefig("milnica.png")

minimalna ploskev

Napolnitev matrike ob eliminaciji

Koda

NumMat.desne_strani — Method.
desne_strani(s, d, z, l)

Izračunaj desne strani pri reševanju robnega problema za Laplaceovo enačbo v 2 dimenzijah.

Argumenti

  • s::Vector: robne vrednosti na spodnjem robu
  • d::Vector: robne vrednosti na desnem robu
  • z::Vector: robne vrednosti na zgornjem robu
  • l::Vector: robne vrednosti na levem robu
source
NumMat.matrika_sistema — Method.
L = matrika_sistema(n, m)

Zapiši matriko za Laplaceov operator v 2D na pravokotnem območju. Matrika L je matrika sistema enačb za diskretizirano laplaceovo enačbo

\[u_{i-1,j}+u_{i,j-1} - 4u_{ij} + u_{i+1,j}+u_{i,j+1} = 0.\]
source
NumMat.resi_robni_problem — Method.
resi_robni_problem(fs, fd, fz, fl, h, meje=((a, b), (c, d)))

Izračunaj približek za rešitev robnega problema za Laplaceovo enačbo

\[\triangle u(x,y) = 0\]

na pravokotniku $[a, b]\times[c, d]$, kjer so vrednosti na robu podane s funkcijami $u(x, c) = f_s(x)$, $u(b, y) = f_d(y)$, $u(x, d) = f_z(x)$ in $u(a, y) = f_l(y)$.

Približek poiščemo z metodo deljenih diferenc s korakom h.

Rezultat

  • Z::Matrix je matrika vrednosti rešitve v notranjosti in na robu.
  • x::Vector je vektor vrednosti abscise
  • y::Vector je vektor vrednosti ordinate

Primer

using Plots
Z, x, y = resi_robni_problem(x->sin(x), y->-sin(y), 
            x->sin(x), y->-sin(x), 0.1, ((0, pi), (0, pi)))
surface(x, y, Z)
source