Eksamen MEK1100 2020¶
Løsningsforslag¶
Dette l?sningsforslaget bruker b?de programmering og direkte regning til ? l?se de fleste oppgavene. Det er strengt tatt ikke n?dvendig ? bruke programmering p? annet enn plotteoppgavene (2b, 2d, 3b) og beregning av buelengde (3d).
Oppgave 1 - Divergens og curl¶
Finn divergensen og virvlingen til f?lgende vektorfelt
- i) $\vec{u} = x \,\mathbf{i} + y \,\mathbf{j} + z \,\mathbf{k}$
- ii) $\vec{u} = r \cos \theta \,\mathbf{i_{r}} + r \sin \theta \,\mathbf{i}_{\theta} + z \,\mathbf{k}$
- iii) $\vec{u} = \mathbf{i_r} + \mathbf{i}$
Her er $\mathbf{i}$, $\mathbf{j}$ og $\mathbf{k}$ de Kartesiske enhetsvektorene i henholdsvis $x$, $y$ og $z$-retning, mens $\mathbf{i}_r, \mathbf{i}_{\theta}, \mathbf{k}$ er enhetsvektorene i sylindriske koordinater for $r, \theta$ og $z$-retning.
import sympy as sp
import numpy as np
from sympy.vector import CoordSys3D
from IPython.display import Math
N = CoordSys3D('N')
x, y, z = sp.symbols('x,y,z', real=True)
ue = (x, y, z)
divu = ue[0].diff(x, 1) + ue[1].diff(y, 1) + ue[2].diff(z, 1)
divu
# Eventuelt
us = N.x*N.i + N.y*N.j + N.z*N.k
sp.vector.divergence(us)
Virvling¶
\begin{align*} \nabla \times \vec{u} &= \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ x & y & z \end{vmatrix} \\ &= \mathbf{i}\left(\frac{\partial z}{\partial y} - \frac{\partial y} {\partial z} \right) + \mathbf{j}\left(\frac{\partial x}{\partial z} - \frac{\partial z} {\partial x} \right) + \mathbf{k}\left(\frac{\partial y}{\partial x} - \frac{\partial x} {\partial y} \right) \\ &= \vec{0} \end{align*}Virvlingen er lik nullvektor.
Implementering¶
curl = (ue[2].diff(y, 1) - ue[1].diff(z, 1), ue[0].diff(z, 1) - ue[2].diff(x, 1), ue[1].diff(x, 1) - ue[0].diff(y, 1))
curl
# Eventuelt
sp.vector.curl(us)
ii)¶
Har n? sylinderkoordinater og
$$ \vec{u} = r \cos \theta \,\mathbf{i_{r}} + r \sin \theta \,\mathbf{i}_{\theta} + z \,\mathbf{k}. $$Divergens¶
Finner divergensen
\begin{align*} \nabla \cdot \vec{u} &= \frac{1}{r}\left(\frac{\partial r u_r}{\partial r} + \frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial r u_z}{\partial z} \right) \\ &= \frac{1}{r}\left(\frac{\partial r^2 \cos \theta}{\partial r} + \frac{\partial r \sin \theta}{\partial \theta} + \frac{\partial r z}{\partial z}\right) \\ &= \frac{1}{r}\left( 2r \cos \theta + r \cos \theta + r\right) \\ &= 3 \cos \theta + 1 \end{align*}Implementering¶
r, theta, z = sp.symbols('x,y,z', real=True, positive=True)
ue = (r*sp.cos(theta), r*sp.sin(theta), z)
divu = 1/r*((ue[0]*r).diff(r, 1) + ue[1].diff(theta, 1) + (ue[2]*r).diff(z, 1))
divu = sp.simplify(divu)
Math('\\text{div}(u) = ' +sp.latex(divu, symbol_names={r: 'r', theta: '\\theta', z: 'z'}))
# Eventuelt
B = N.create_new('B', transformation='cylindrical')
us = B.r*sp.cos(B.theta)*B.i + B.r*sp.sin(B.theta)*B.j + B.z*B.k
divu = sp.vector.divergence(us)
Math('\\text{div}(u) = ' + sp.latex(divu).replace('{theta}_{B}', '\\theta'))
Virvling¶
\begin{align*} \nabla \times \vec{u} &= \frac{1}{r} \begin{vmatrix} \mathbf{i}_r & r\mathbf{i}_{\theta} & \mathbf{i}_z \\ \frac{\partial }{\partial r} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ r \cos \theta & r^2 \sin \theta & z \end{vmatrix} \\ &= \frac{1}{r}\left( \mathbf{i}_r \left( \frac{\partial z}{\partial \theta} - \frac{\partial r^2 \sin \theta}{\partial z} \right) +r \mathbf{i}_{\theta} \left( \frac{\partial r \cos \theta}{\partial z} - \frac{\partial z}{\partial r}\right) +\mathbf{i}_z \left( \frac{\partial r^2 \sin \theta}{\partial r} - \frac{\partial r \cos \theta}{\partial \theta} \right) \right) \\ &= \frac{1}{r}\mathbf{i}_z\left(2 r \sin \theta + r \sin \theta \right) \\ &= 3 \sin \theta \, \mathbf{i}_z \end{align*}Implementering¶
curl = (1/r*(ue[2].diff(theta, 1) - (r*ue[1]).diff(z, 1)),
(ue[0].diff(z, 1) - ue[2].diff(r, 1)),
(1/r)*((r*ue[1]).diff(r, 1) - ue[0].diff(theta, 1)))
curl = sp.simplify(curl)
Math('\\text{curl}(u) = ' + sp.latex(curl, symbol_names={theta: '\\theta'}))
# Eventuelt
curlu = sp.vector.curl(us)
Math('\\text{curl}(u) = ' + sp.latex(curlu).replace('{theta}_{B}', '\\theta').replace('\\hat{k}_{B}', '\\mathbf{k}' ))
iii)¶
$$ \vec{u} = \mathbf{i_r} + \mathbf{i} $$Det er flere m?ter ? l?se denne oppgaven p?. Man kan velge ? transformere vektoren til Kartesiske eller sylinder-koordinater, men det er ikke n?dvendig da b?de divergens og virvling er line?re operatorer.
Ved konvertering kan man bruke
$$ \mathbf{i}_r = \cos \theta \,\mathbf{i} + \sin \theta \,\mathbf{j}, $$og
$$ r = \sqrt{x^2+y^2}, \quad x = r \cos \theta, \quad y = r \sin \theta, $$og skrive vektoren i Kartesiske koordinater som
$$ \vec{u} = \left(1+\frac{x}{\sqrt{x^2+y^2}} \right) \mathbf{i} + \frac{y}{\sqrt{x^2+y^2}} \,\mathbf{j}. $$Alternativt, bruk
$$ \mathbf{i} = \cos \theta \mathbf{i}_r - \sin \theta \mathbf{i}_{\theta}, $$og skriv vektoren som
$$ \vec{u} = \left( 1 + \cos \theta \right) \mathbf{i}_r - \sin \theta \mathbf{i}_{\theta}. $$Alle fremgangsm?ter skal gi samme svar.
Divergens¶
Eventuelt
\begin{align*} \nabla \cdot (\mathbf{i}_r + \mathbf{i}) &= \nabla \cdot \mathbf{i}_r + \nabla \cdot \mathbf{i}\\ &= \nabla \cdot \mathbf{i}_r = \frac{1}{r} \frac{\partial r}{\partial r} = \frac{1}{r} \end{align*}Eller (hopper over litt lang, men triviell utledning.)
\begin{align*} \nabla \cdot \vec{u} &= \frac{\partial }{\partial x}\left(1+\frac{x}{\sqrt{x^2+y^2}} \right) + \frac{\partial}{\partial y}\left(\frac{y}{\sqrt{x^2+y^2}}\right) \\ &= \ldots \\ &= \frac{1}{\sqrt{x^2+y^2}}\\ &= \frac{1}{r} \end{align*}Til slutt ved rene sylinderkoordinater
\begin{align*} \nabla \cdot \vec{u} &= \frac{1}{r} \frac{\partial r(1+\cos \theta)}{\partial r} - \frac{1}{r}\frac{\partial \sin \theta}{\partial \theta} \\ &= \frac{1+\cos \theta}{r} - \frac{\cos \theta}{r} \\ &= \frac{1}{r} \end{align*}Implementering¶
ue = (1+x/sp.sqrt(x**2+y**2), y/sp.sqrt(x**2+y**2), sp.S.Zero)
divu = ue[0].diff(x, 1) + ue[1].diff(y, 1) + ue[2].diff(z, 1)
divu = sp.simplify(divu)
divu
# Eventuelt
us = (1+N.x/sp.sqrt(N.x**2+N.y**2))*N.i + N.y/sp.sqrt(N.x**2+N.y**2)*N.j
divu = sp.vector.divergence(us)
divu = sp.simplify(divu)
divu
Her er $x_N$ og $y_N$ coordinatene i coordinatsystem N
, alts? $x$ og $y$. Kunne ogs? ha implementert i sylinder-koordinater.
Virvling¶
Eventuelt
\begin{align*} \nabla \times (\mathbf{i}_r + \mathbf{i}) &= \nabla \times \mathbf{i}_r + \nabla \times \mathbf{i}\\ &= \nabla \times \mathbf{i}_r = \vec{0} \end{align*}siden kun konstant koeffisient og $h_r=1$.
Eller
\begin{align*} \nabla \times \vec{u} &= \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ 1+\frac{x}{r} & \frac{y}{r} & 0 \end{vmatrix} \\ &= \mathbf{i}\left(\frac{\partial 0}{\partial y} - \frac{\partial y/r} {\partial z} \right) + \mathbf{j}\left(\frac{\partial 1+x/r}{\partial z} - \frac{\partial 0} {\partial x} \right) + \mathbf{k}\left(\frac{\partial y/r}{\partial x} - \frac{\partial 1+x/r} {\partial y} \right) \\ &= \mathbf{k} \left( -\frac{xy}{r^3} + \frac{xy}{r^3} \right) \\ &= \vec{0} \end{align*}Eller
\begin{align*} \nabla \times \vec{u} &= \frac{1}{r}\begin{vmatrix} \mathbf{i}_r & r\mathbf{i}_{\theta} & \mathbf{k} \\ \frac{\partial }{\partial r} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z} \\ 1+\cos \theta & \sin \theta & 0 \end{vmatrix} \\ &= \frac{1}{r} \mathbf{i}_z \left(-\sin \theta - (-\sin \theta) \right)\\ &= \vec{0} \end{align*}Implementering¶
curl = (ue[2].diff(y, 1) - ue[1].diff(z, 1), ue[0].diff(z, 1) - ue[2].diff(x, 1), ue[1].diff(x, 1) - ue[0].diff(y, 1))
curl
# Eventuelt
sp.vector.curl(us)
Oppgave 2 - Elliptiske koordinater¶
Elliptiske koordinater $(u, v)$ er gitt ved posisjonsvektor
$$ \vec{r} = a \cosh u \cos v \mathbf{i} + a \sinh u \sin v \mathbf{j}. $$2a) Finn enhetsvektorer og skaleringsfaktorer.¶
- Finn enhetsvektorene $\mathbf{e}_{u}$, $\mathbf{e}_{v}$ og skaleringsfaktorene til de elliptiske koordinatene. Er enhetsvektorene ortogonale?
Skaleringsfaktorene er gitt ved
$$ h_u = \left|\frac{\partial \vec{r}}{\partial u}\right|, \quad h_v = \left|\frac{\partial \vec{r}}{\partial v}\right|. $$Enhetsvektorene er gitt ved
$$ \mathbf{e}_u = \frac{1}{h_u}\frac{\partial \vec{r}}{\partial u}, \quad \mathbf{e}_v = \frac{1}{h_v}\frac{\partial \vec{r}}{\partial v}.$$Finner
\begin{align*} \frac{\partial \vec{r}}{\partial u} &= \frac{\partial \cosh u \cos v}{\partial u} \mathbf{i} + \frac{\partial \sinh u \sin v}{\partial u} \mathbf{j} \\ &= \sinh u \cos v \mathbf{i} + \cosh u \sin v \mathbf{j} \end{align*}og
\begin{align*} \frac{\partial \vec{r}}{\partial v} &= \frac{\partial \cosh u \cos v}{\partial v} \mathbf{i} + \frac{\partial \sinh u \sin v}{\partial v} \mathbf{j} \\ &= -\cosh u \sin v \mathbf{i} + \sinh u \cos v \mathbf{j} \end{align*}hvilket gir, hvis vi bruker $\cosh^2 u = 1 + \sinh^2 u$,
\begin{align*} h_u^2 &= (\sinh u \cos v \mathbf{i} + \cosh u \sin v \mathbf{j}) \cdot (\sinh u \cos v \mathbf{i} + \cosh u \sin v \mathbf{j}) \\ &= \sinh^2 u \cos^2 v + \cosh^2 u \sin^2 v \\ &= \sinh^2 u \cos^2 v + (1+\sinh^2 u) \sin^2 v \\ &= \sinh^2 u (\cos^2 v + \sin^2 v) + \sin^2 v \\ &= \sinh^2 u + \sin^2 v \end{align*}Eventuelt kan vi bruke $\sinh^2 u = \cosh^2 u - 1$ og $\sin^2 v + \cos^2 v = 1$ og f?
\begin{align*} h_u^2 &= \sinh^2 u + \sin^2 v\\ &= \cosh^2 u - 1 + 1 - \cos^2 v \\ &= \cosh^2 u - \cos^2 v \end{align*}Samme svar f?r vi for $h_v$
\begin{align*} h_v^2 &= (-\cosh u \sin v \mathbf{i} + \sinh u \cos v \mathbf{j}) \cdot (-\cosh u \sin v \mathbf{i} + \sinh u \cos v \mathbf{j}) \\ &= \cosh^2 u \sin^2 v + \sinh^2 u \cos^2 v \\ &= \sinh^2 u + \sin^2 v \\ &= \cosh^2 u - \cos^2 v \end{align*}Man kan ogs? bruke andre forenklinger. For eksempel
\begin{equation*} \sinh^{2}u = \frac{1}{2}\left(\cosh 2u - 1 \right), \end{equation*}og da f?r man
\begin{equation*} h_u = h_v = \sqrt{2 \cosh 2u - 2 \cos 2v}. \end{equation*}S? skaleringsfaktorene er
\begin{equation*} h = h_u = h_v = \sqrt{\sinh^2 u \cos^2 v + \cosh^2 u \sin^2 v}. \end{equation*}Med forenkling:
\begin{equation*} h = h_u = h_v = \sqrt{\sinh^2 u + \sin^2 v} = \sqrt{\cosh^2 u - \cos^2 v} = \sqrt{2 \cosh 2u - 2 \cos 2v}. \end{equation*}Alle disse svarene er like gyldige.
Enhetsvektorene er dermed
\begin{align*} \mathbf{e}_u &= \frac{1}{h}\left(\sinh u \cos v \,\mathbf{i} + \cosh u \sin v \,\mathbf{j}\right), \\ \mathbf{e}_v &= \frac{1}{h}\left(-\cosh u \sin v \,\mathbf{i} + \sinh u \cos v \,\mathbf{j}\right). \end{align*}Implementering¶
Starter med ? importere funksjonalitet fra sympy
, og lager to tupler (Python immutable list) for psi=(u, v)
og rv=(sinh(u)*cos(v), sinh(u)*sin(v))
import sympy as sp
import numpy as np
u, v = psi = sp.symbols('u,v', real=True)
rv = (sp.cosh(u)*sp.cos(v), sp.sinh(u)*sp.sin(v))
Bruker noen mer eller mindre generelle funksjoner for a finne basisvektorer, skaleringsfaktorer og enhetsvektorer. Disse er mer eller mindre eksakt som allerede g?tt igjennom p? forelesning for parabolske koordinater:
def basisvektorer(psi, rv):
"""Returner basisvektorer
Parameters
----------
psi : Tuple av nye variable
rv : Posisjonsvektor
"""
b = np.zeros((len(psi), len(rv)), dtype=object)
for i, ui in enumerate(psi):
for j, rj in enumerate(rv):
b[i, j] = sp.simplify(rj.diff(ui, 1))
return b
def skaleringsfaktorer(b):
"""Returner skaleringsfaktorer
Parameters
----------
b : basisvektorer
"""
h = np.zeros(b.shape[0], dtype=object)
for i, s in enumerate(np.sum(b**2, axis=1)):
h[i] = sp.simplify(sp.sqrt(s))
# Sympy er dessverre ikke smart nok til ? forenkle videre. Hjelper til med denne:
h[i] = sp.simplify(h[i].replace(sp.cosh(u)**2, 1+sp.sinh(u)**2))
return h
def enhetsvektorer(psi, rv):
"""Returner enhetsvektorer og skaleringsfaktorer
Parameters
----------
psi : Tuple av nye variable
rv : Posisjonsvektor
"""
b = basisvektorer(psi, rv)
hi = skaleringsfaktorer(b)
return b / hi[None, :], hi
e, hi = enhetsvektorer(psi, rv)
Merk at vi velger l?sningen med $h = \sqrt{\sinh^2 u + \sin^2 v}$.
hi = sp.sympify(hi)
hi
Skriver ut enhetsvektorene etter f?rst ? forenkle noe s? det ikke blir s? rotete:
h = sp.Symbol('h')
replace = lambda f: f.replace((sp.sin(v)**2+sp.sinh(u)**2)**(-1/2), 1/h)
es = replace(sp.sympify(e))
es
Printer ut sammen med de Kartesiske enhetsvektorene
Math('\\mathbf{e_u} = ' + sp.latex(es[0, 0]) + '\\mathbf{i} +' + sp.latex(es[0, 1]) + '\\mathbf{j} \\\\'
+'\\mathbf{e_v} = ' + sp.latex(es[1, 0]) + '\\mathbf{i} +' + sp.latex(es[1, 1]) + '\\mathbf{j}')
Ortogonalitet¶
Enhetsvektorene er ortogonale siden
\begin{equation*} \mathbf{e}_u \cdot \mathbf{e}_v = 0, \end{equation*}og
\begin{equation*} \mathbf{e}_u \cdot \mathbf{e}_u = \mathbf{e}_v \cdot \mathbf{e}_v = 1. \end{equation*}Har i koden at $e[0]=\mathbf{e}_u$ og $e[1]=\mathbf{e}_v$. Bruker dette:
np.sum(e[0]*e[1])
sp.simplify(np.sum(e[0]*e[0]))
Her ser vi at sympy trenger litt hjelp til ? forenkle
sp.simplify((np.sum(e[0]*e[0])).replace(sp.cosh(u)**2, 1+sp.sinh(u)**2))
Den er med andre ord ortogonal. Samme med $\mathbf{e}_v$:
sp.simplify((np.sum(e[1]*e[1])).replace(sp.cosh(u)**2, 1+sp.sinh(u)**2))
2b) Operatorer¶
- Gitt et skalarfelt $f$ og en vektor $\vec{w} = w_u \mathbf{e}_{u} + w_v \mathbf{e}_v$. Gi $\nabla f$, $\nabla \cdot \vec{w}$ og Laplace operatoren $\nabla^2$ i elliptiske koordinater.
Under bruker vi $h = \sqrt{\sinh^2 u + \sin^2 v } = \sqrt{\cosh^2 u - \cos^2 v} = \sqrt{2 \cosh 2u - 2 \cos 2v}$ i formler fra Vector Calculus. Merk at man m? forenkle for $h_u=h_v$ for ? f? full pott.
$$ \nabla f = \frac{\mathbf{e}_u}{h}\frac{\partial f}{\partial u} + \frac{\mathbf{e}_v}{h}\frac{\partial f}{\partial v}, $$$$ \nabla \cdot \vec{w} = \frac{1}{h^2}\left(\frac{\partial w_u h}{\partial u} + \frac{\partial w_v h}{\partial v} \right), $$$$ \nabla^2 = \frac{1}{h^2}\left( \frac{\partial^2}{\partial u^2} + \frac{\partial^2}{\partial v^2}\right). $$For divergensen kan man ekspandere de deriverte, men uttrykket blir ikke enklere.
2c) Skisser koordinatkurvene¶
- Presenter en skisse av koordinatkurvene i et Kartesisk koordinatsystem. Det vil si, skisser kurver med konstante verdier for b?de $u$ eller $v$, jamf?r Fig. 6.3 i Vector Calculus.
import matplotlib.pyplot as plt
%matplotlib inline
M = 32
# Velger 3 ganger fler punkter i v retning for ? f? jevne intervall
ui = np.broadcast_to(np.linspace(0, 1, M)[:, None], (M, 3*M))
vi = np.broadcast_to(np.linspace(0, 2*np.pi, 3*M)[None, :], (M, 3*M))
# Varierer c mellom 0 og 1 og veksler mellom ? holde u og v konstant
for c in np.linspace(0, 1, 10):
plt.plot(np.cosh(c)*np.cos(vi[0]), np.sinh(c)*np.sin(vi[0]), 'b') # Her er c en konstant u
plt.plot(np.cosh(ui[:, 0])*np.cos(np.pi*c), np.sinh(ui[:, 0])*np.sin(np.pi*c), 'r') # Konstant v (?verste halvdel - 0 til pi)
plt.plot(np.cosh(ui[:, 0])*np.cos(np.pi*(c+1)), np.sinh(ui[:, 0])*np.sin(np.pi*(c+1)), 'r') # Konstant v (nederste halvdel - pi til 2pi)
plt.xlabel('x')
plt.ylabel('y')
2d) Kontur- og pilplott¶
Lag konturplott av skalarfeltet $$ f(u, v) = (1-u^2)\cos(2v), $$ i b?de elliptiske og Kartesiske koordinater.
Finn $\nabla f$ (i elliptiske eller Kartesiske koordinater) og lag et pilplott av denne i samme figur som det Kartesiske konturplottet. Kommenter retningen p? pilene.
Lager skalarfelt $f(u, v) = (1-u^2) \cos(2v)$
f = (1-u**2)*sp.cos(2*v)
Merk at vi bruker $x=2uv$ og $y=u^2-v^2$ evaluert p? et strukturert grid. sp.lambdify
er en effektiv (vektorisert) metode ? evaluere en sympy
funksjon p?. S? under tilsvarer f(u, v) = sp.lambdify((u, v), f)(ui, vi)
.
fj = sp.lambdify((u, v), f)(ui, vi)
Hvis vi n? velger ? plotte $f(u, v)$ i det nye koordinatsystemet f?r vi.
plt.contourf(ui, vi, fj)
plt.xlabel('u')
plt.ylabel('v')
Men det er kanskje mer interessant ? se resultatet i fysiske (Kartesiske) koordinater. Vi trenger derfor ? finne kartesiske x, y
fra de gitte u, v
. Gj?r dette som f?lger
mesh = []
for rj in rv:
mesh.append(sp.lambdify((u, v), rj)(ui, vi))
x, y = mesh
plt.contourf(x, y, fj)
plt.xlabel('x')
plt.ylabel('y')
? plotte gradienten i Kartesiske koordinater er mer involvert siden vi har beregnet gradienten i de nye koordinatene og derfor trenger ? projisere ved ? ta prikk-produktet av gradientvektoren
\begin{align*} \frac{\partial f}{\partial x} &= \nabla f \cdot \mathbf{i},\\ \frac{\partial f}{\partial y} &= \nabla f \cdot \mathbf{j}. \end{align*}For ? finne gradientvektoren deriverer vi f?rst for ? finne komponentene til $\nabla f$ i nye koordinater
df = np.array((1/hi[0]*f.diff(u, 1), 1/hi[1]*f.diff(v, 1)))
Merk at df
n? ikke inneholder enhetsvektorer. S? f?r vi prikker med $\mathbf{i}$ og $\mathbf{j}$ m? vi gange med enhetsvektorene $\mathbf{e_u}$ og $\mathbf{e_v}$ for ? f? $\nabla f$
Math('\\nabla f = ' + sp.latex(replace(df[0])) + '\mathbf{e}_u' + sp.latex(replace(df[1])) + '\mathbf{e}_v')
Det kan v?re morsomt ? se p? denne ogs? i Kartesiske koordinater selv om vi bare sp?r etter $\nabla f$ i elliptiske eller Kartesiske koordinater. Da m? vi multiplisere med enhetsvektorene
gradf = e[0]*df[0] + e[1]*df[1]
gradf = sp.sympify(gradf)
Merk at vi med denne summen n? har f?tt satt inn for $\mathbf{e_u}$ og $\mathbf{e_v}$, s? vektoren gradf
over er gitt ved Kartesiske enhetsvektorer.
Math('\\nabla f = ' + sp.latex(sp.simplify(gradf[0]))+'\mathbf{i} \\\\ ' + sp.latex(sp.simplify(gradf[1]))+'\mathbf{j}')
Ved prikking av denne vektoren mot $\mathbf{i}$ er resultatet gradf[0]
, mens prikking mot $\mathbf{j}$ gir gradf[1]
. Derfor hopper vi over direkte regning av prikkproduktet og henter ganske enkelt de Kartesiske vektorkomponentene vi beh?ver til plottet
import warnings
warnings.filterwarnings("ignore")
dfdxi = sp.lambdify((u, v), gradf[0])(ui, vi)
dfdyi = sp.lambdify((u, v), gradf[1])(ui, vi)
plt.contourf(x, y, fj)
plt.quiver(x[::2, ::2], y[::2, ::2], dfdxi[::2, ::2], dfdyi[::2, ::2], pivot='mid')
plt.xlabel('x')
plt.ylabel('y')
Merk at gradienten peker i retning av ?kende $f$.
Oppgave 3¶
En Taylor-Green virvel er en to-dimensjonell analytisk modell for periodiske str?mningsvirvler som avtar i styrke over tid. Virvlene er gitt ved
$$ \vec{u}(x,y,t) = \left(\cos x \sin y \,\mathbf{i} - \sin x \cos y \,\mathbf{j} \right) \exp^{-2\nu t}, $$for tiden $t>=0$ i et domene $\Omega = [-\pi, \pi] \times [-\pi, \pi]$. Parameteren $\nu$ er ? regne som en konstant.
3a) Strømfunksjon og skalarpotensial¶
Vis at str?mfunksjonen er
$$ \psi(x, y, t) = \cos x \cos y \exp^{-2 \nu t}. $$
Hvordan kunne vi vite p? forh?nd at dette feltet har en str?mfunksjon?
Vi sjekker at divergensen til $\vec{u}$ er lik 0.
\begin{align*} \nabla \cdot \vec{u} &= \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \\ &= \left(-\sin x \sin y - \sin x (-\sin y) \right) \exp^{-2\nu t} \\ &= 0 \end{align*}Vi finner $\psi$ enkelt ved ? l?se disse
\begin{align*} -\frac{\partial \psi}{\partial y} &= \cos x \sin y \exp^{-2\nu t} \\ \frac{\partial \psi}{\partial x} &= -\sin x \cos y \exp^{-2\nu t} \end{align*}Integrerer f?rste og andre likning mhp $y$ og $x$ og f?r
\begin{align*} \psi &= \cos x \cos y \exp^{-2\nu t} + f_1(x) \\ \psi &= \cos x \cos y \exp^{-2\nu t} + f_2(y) \\ \end{align*}Ser at $f_1(x) = f_2(y) = $ konstant. Vi kan velge konstanten fritt for en str?mfunksjon og velger 0. Da har vi vist at str?mfunksjonen er som gitt i starten av oppgaven.
Implementering¶
nu, t = sp.symbols('nu t', real=True, positive=True)
us = (sp.cos(N.x)*sp.sin(N.y)*N.i - sp.sin(N.x)*sp.cos(N.y)*N.j)*sp.exp(-2*nu*t)
divu = sp.vector.divergence(us)
divu
- Har dette vektorfeltet et skalarpotensial? Hvis ja, finn skalarpotensialet.
Vi sjekker om vektorfeltet er virvelfritt. Vi har et 2D vektorfelt, s? virvlingen kan kun ha én komponent, i $z$-retning
\begin{align*} \nabla \times \vec{v} &= \mathbf{k} \left(\frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \right)\\ &= \mathbf{k} \left( -\cos x \cos y - \cos x \cos y \right) \exp^{-2 \nu t} \\ &= -2 \mathbf{k} \cos x \cos y \exp^{-2 \nu t} \end{align*}Implementering¶
curlu = sp.vector.curl(us)
curlu
Siden feltet ikke er virvelfritt har det ikke et skalarpotensial.
3b) Pilplott med strømlinjer¶
- Lag et pilplott av vektorfeltet $\vec{u}(x, y, 0)$. Plott i samme figur str?mlinjene til $\vec{u}$.
M = 30
xi, yi = np.meshgrid(*(np.linspace(-np.pi, np.pi, M),)*2)
uss = us.subs({t: 0, nu: 1})
usx, usy = uss.dot(N.i), uss.dot(N.j)
ux = sp.lambdify((N.x, N.y), usx)(xi, yi)
uy = sp.lambdify((N.x, N.y), usy)(xi, yi)
plt.quiver(xi, yi, ux, uy)
psi = sp.cos(N.x)*sp.cos(N.y)
psij = sp.lambdify((N.x, N.y), psi)(xi, yi)
plt.contour(xi, yi, psij)
plt.xlabel('x')
plt.ylabel('y')
3c) Fluks og sirkulasjon¶
Hva blir fluksen
$$ \oint_{C} \vec{u} \cdot \vec{n} ds, $$
av vektorfeltet $\vec{u}$ ut av et rektangul?rt omr?de $\Omega = [0, \pi/2] \times [0, \pi/2]$ som omsluttes av kurven $C$? Her er $\vec{n}$ normalvektoren som peker ut fra omr?det og $ds$ er et linjeelement langs kurven $C$. Forklar hvorfor man kan bruke Gauss divergensteorem her selv om det bare er et to-dimensjonalt integral.
Merk at Gauss divergensteorem normalt gjelder en sammenheng mellom et flateintegral og et volumintegral.
Man kan vise at Gauss teorem gjelder i 2D ved ? vise til at det er en generalisering av Greens teorem, som forklart i kap 6.14.3 i Linstr?m og Hveberg "Flervariabel analyse med line?r algebra".
Man kan ogs? forestille seg at 2D kvadratet utvides til en kube ved ? ekstrahere i $z$-retningen en lengde $L$ (lengden er ikke viktig, bare at domenet n? blir 3D). Fluksintegralet over denne kuben er
$$ \oint_{A} \vec{u} \cdot \vec{n} \,d\sigma, $$der $d\sigma$ er et flateelement og $A$ er den omsluttende flaten til kuben. Vi vet at siden $\vec{u}$ ikke er en funksjon av $z$ s? vil, for dette fluks-integralet, de to flateintegralene i $xy$-planet kansellere (flatene ved $z=0$ og $z=L$). S? vi st?r igjen med et flateintegral rundt kurven $C$, ekstrahert en dybde L.
$$ \oint_{A} \vec{u} \cdot \vec{n} d\sigma = L \oint_{C} \vec{u} \cdot \vec{n} \,ds. $$Venstresiden kan vi bruke vanlig Gauss teorem p?, s?
$$ \int_V \nabla \cdot \vec{u} \, dV = L \oint_{C} \vec{u} \cdot \vec{n} \,ds. $$Dermed er det klart at vi kan bruke vanlig Gauss teorem til ? finne integralet p? h?yresiden. Vi kan eventuelt forenkle videre med ? bruke $dV = dx dx dz$
\begin{align*} \int_{0}^L \int_{0}^{\pi/2} \int_0^{\pi/2} \nabla \cdot \vec{u} \, dxdydz &= \int_{0}^{L} dz \int_0^{\pi/2} \int_0^{\pi/2} \nabla \cdot \vec{u} \, dxdy, \\ &= L \int_0^{\pi/2} \int_0^{\pi/2} \nabla \cdot \vec{u} \, dxdy, \end{align*}s? vi blir kvitt $L$'en. Men dette er ikke n?dvendig. Vi vet allerede at $\nabla \cdot \vec{u}=0$, s? vi ender opp med at
$$ \oint_{C} \vec{u} \cdot \vec{n} ds = 0. $$Hva blir sirkulasjonen
$$ \oint_{C} \vec{u} \cdot d\vec{r}, $$ n?r vi beveger oss mot klokka (alts? fra $(0, 0)$ til $(\pi/2, 0)$ osv.)? Finn resultatet b?de ved direkte regning av kurveintegralet, og ved ? benytte et passende integralteorem.
Ved direkte regning deler vi opp linjeintegralet over de fire siden av kvadratet:
\begin{align*} \oint_{C} \vec{u} \cdot d\vec{r} &= \int_{x=0}^{\pi/2} \vec{u}(x, 0) \cdot \mathbf{i} \,dx \\ &+ \int_{y=0}^{\pi/2} \vec{u}(\pi /2, y) \cdot \mathbf{j} \,dy \\ &+ \int_{x=\pi/2}^{0} \vec{u}(x, \pi/2) \cdot \mathbf{i} \,dx\\ &+ \int_{y=\pi/2}^{0} \vec{u}(0, y) \cdot \mathbf{j} \,dy. \end{align*}Merk at vi integrerer motsatt vei enten ved ? prikke mot negativ enhetsvektor, eller ved ? snu integrasjonsgrensene, men ikke begge to. Da kanselleres effekten. Integralet er rett frem og vi kan bruke sympy:
(sp.integrate(us.dot(N.i).subs(N.y, 0), (N.x, 0, sp.pi/2))+
sp.integrate(us.dot(N.j).subs(N.x, sp.pi/2), (N.y, 0, sp.pi/2))+
sp.integrate(us.dot(N.i).subs(N.y, sp.pi/2), (N.x, sp.pi/2, 0))+
sp.integrate(us.dot(N.j).subs(N.x, 0), (N.y, sp.pi/2, 0)))
Ved ? benytte Stokes' teorem kan vi gjenbruke curlen:
$$ \oint_{C} \vec{u} \cdot d\vec{r} = \int_A (\nabla \times \vec{u}) \cdot \mathbf{k} \, d\sigma. $$sp.integrate(sp.integrate(curlu.dot(N.k), (N.x, 0, sp.pi/2)), (N.y, 0, sp.pi/2))
Siden dette vektorfeltet er i 2D, s? kan man ogs? velge ? bruke Green's teorem, siden Stokes er en generalisering av Green's teorem til 3 dimensjoner.
3d) Ekviskalarflate og buelengde¶
- La $\psi(x, y, t=0)=z$ representere h?yde og $\beta(x, y, z) = z - \cos x \cos y = 0$ en ekviskalarflate. Finn flatenormalen.
Normaliserer (see s 28-29 i Gjevik) slik at flatenormalen $\vec{n}$ er gitt ved
$$ \vec{n} = \frac{\nabla \beta}{|\nabla \beta|}, $$$$ \vec{n} = \frac{\sin x \cos y \mathbf{i} + \cos x \sin y \mathbf{j} + \mathbf{k}}{\sqrt{\sin^2 x \cos^2 y + \cos^2 x \sin^2 y + 1}}. $$Her kan nevner forenkles videre med ? bruke $\cos^2 x = 1-\sin^2 x$, men det er ikke n?dvendig. Sympy gj?r denne 'forenklingen':
beta = N.z - psi
n = sp.vector.gradient(beta)
n = sp.simplify(n / n.magnitude())
n
- Hvis man holder seg p? ekviskalarflaten $\beta$ og g?r en full sirkel med radius 1 ($x^2+y^2=1$) rundt origo, hvor lang er da denne buelengden?
Denne oppgaven er veldig utmattende ? l?se med penn og papir, men kan l?ses enkelt i Geogebra ved ? bruke samme fremgangsm?te som i denne appleten, som ble g?tt igjennom i forelesning. L?ses kanskje enda enklere i sympy. Har parametrisering
$$ \vec{r} = \cos(\theta)\mathbf{i} + \sin(\theta) \mathbf{j} + \cos(\cos(\theta)) \cos(\sin(\theta)) \mathbf{k}, $$og buelengden er
$$ \int_0^{2\pi} \left|\frac{\partial \vec{r}}{\partial \theta}\right| d\theta. $$I sympy:
theta = sp.Symbol('theta', real=True, positive=True)
r = sp.cos(theta)*N.i + sp.sin(theta)*N.j + sp.cos(sp.cos(theta))*sp.cos(sp.sin(theta))*N.k
dr = r.diff(theta, 1)
D = sp.Integral(dr.magnitude(), (theta, 0, 2*sp.pi)).evalf(20) # m? evaluere numerisk
print('Buelengden =', D)
Vi beregnet buelengde i Python og Sympy i dette scriptet som ble g?tt igjennom i forelesning.
3e) Finn trykket¶
Newton's andre lov for inkompressibel str?mning (inkludert friksjon) gir oss momentumlikningen
$$ \frac{D \vec{u}}{D t} = -\nabla p + \nu \nabla^2 \vec{u}. \label{eq:NS} $$
Vis at den partikkelderiverte av feltet $\vec{u}$ er gitt ved
dudt = us.diff(t, 1)
dudt
conv = sp.simplify(us.dot(sp.vector.Del())(us))
conv
Ser at dudt
$=-2\nu \vec{u}$ og conv
$=- 1/2 \left( \sin 2 x \mathbf{i} + \sin 2 y \mathbf{j} \right) \exp^{-4 \nu t}$.
Bruk dette og videre innsetting i momentumlikningen til ? finne trykket $p$.
Finner f?rst Laplace uttrykket:
nu*sp.vector.Laplacian(us).doit()
Ser at Laplace og tidsderivert kansellerer eksakt:
sp.simplify(nu*sp.vector.Laplacian(us).doit() - dudt)
St?r igjen med
$$ (\vec{u} \cdot \nabla) \vec{u} = -\nabla p. $$Har da
\begin{align*} \frac{\partial p}{\partial x} &= -(\vec{u} \cdot \nabla) \vec{u} \cdot \mathbf{i}, \\ \frac{\partial p}{\partial y} &= -(\vec{u} \cdot \nabla) \vec{u} \cdot \mathbf{j}. \end{align*}Eller
\begin{align*} \frac{\partial p}{\partial x} &= -\frac{1}{2}\sin 2x \exp^{-4 \nu t}, \\ \frac{\partial p}{\partial y} &= -\frac{1}{2}\sin 2y \exp^{-4 \nu t}. \end{align*}Som kan l?ses til
$$ p(x, y, t) = -\frac{1}{4}(\cos 2x + \cos 2y) \exp^{-4 \nu t} + p_0, $$der $p_0$ er en vilk?rlig konstant.