IN-KJM1900 - H?st 2018

Koblede diffligninger og interpolasjon

IN-KJM1900, Fjerde forelesning, 6.november

Mål for økta:

  • Oversikt over pensum
  • Introduksjon av deloppgave 2
  • Detaljert gjennomgang av hydrologisk delmodell
  • Ny teori:
    • Koblede differensialligninger
    • Interpolasjon

Merk: Vi rakk ikke ? gjennomg? hele den hydrologiske delmodellen, s? deler av dette vil bli gjennomg?tt 13.november.

Orakeltjeneste, uke 45

Orakelhjelp tilbys ved Kjemisk Institutt p? f?lgende tidspunkter denne uka

Tirsdag 6/11 14.00-16.00
Onsdag 7/11 13.00-15.00
Torsdag 8/11 13.00-15.00

Merk deg at torsdags?kta er flyttet, og at den ikke lenger er i konflikt med samretting fra 10.15-12.00.

Se ogs? nettressurser.

Innlevering av deloppgave 2

  • Fristen for ? levere deloppgave 2 av prosjektet er s?ndag 18. november klokka 23.59.
  • Oppgaven leveres (som normalt) p? Devilry.?
    • ?
  • Vi oppfordrer til ? benytte samrettinga p? torsdagene klokka 10.15-12.00
  • Gruppetimene g?r som normalt
  • Sp?rsm?l om prosjektet?

Pensum IN-KJM1900

(Vil ogs? publiseres p? semestersidene.)

Felles med IN1900, fra A Primer on Scientific Programming with Python, Fifth Edition av Hans Petter Langtangen:

Kapittel Innhold Unntak
1 Computing with formulas 1.6 og 1.7
2 Loops and lists Ingen
3 Functions and branching 3.3
4 User input and error handling 4.4, 4.8, 4.10
5 Array computing and curve plotting 5.5.1, 5.5.3, 5.7 og 5.9-5.12

Se for?vrig ogs? semestersidene til IN1900.

Pensum IN-KJM1900

Fra samme bok, men dekket i kjemidelen av kurset:

Kapittel Innhold Unntak
6 Dictionaries and strings 6.3-6.7
Appendiks A Sequences and difference equations A.2
Appendiks B Discrete Calculus B.4
Appendiks C Differential equations C.3-C.4
Appendiks E Programming of differential equations E.1.7, E.2.4, E.3-E.4

I tillegg er prosjektoppgaven pensum. Eksamen vil i stor grad sammenfalle med det som introduseres i prosjektet.

Diskretisering

Vi har tidligere sett p? diskret kalkulus, alts? en form for matematikk som beskriver diskrete funksjoner:
Kontinuerlig Diskret
$x$ $x_i = i\cdot\Delta x$
$f(x)$ $f_i = f(x_i) = f(i\cdot\Delta x)$
$\frac{d}{dx}$ $D^+$
$\int$ $\sum$

N?r en funksjon eller variabel er diskretisert, har vi mistet informasjon om verdiene mellom den valgte oppl?sningen. Kan vi finne de igjen?

Interpolasjon

Interpolasjon er en teknikk hvor vi benytter diskrete datapunkter for ? gj?re et mer eller mindre kvalifisert gjett p? hva som finnes mellom.

Dette er p? en m?te en "anti-diskretisering".

Interpolasjon

I dette kurset vil vi benytte en interpolasjons-funksjon fra scipy, som importeres og benyttes p? f?lgende m?te:

In [5]:
from scipy.interpolate import interp1d 
import numpy as np

# data punkter for testing
data_x = np.array([1,2,3,4])
data_y = np.array([1,7,3,5])

y_interp = interp1d(data_x, data_y, 2) #gj?r en andre ordens interpolasjon

x = np.linspace(1,4,100) #?k oppl?sningen p? grid
y = y_interp(x) #beregn funksjonsverdiene til den interpolerte funksjonen

Kopier gjerne koden over og kj?r p? din egen maskin for ? pr?ve den ut.

Diskusjon: hvor sikre kan vi v?re p? de interpolerte verdiene mellom datapunktene?

$\rightarrow$ Kodeeksempel aliasering

Prosjektets del 2: Birkenesmodellen

  • Numerisk modell for vannbevegelse og kjemi i jordv?ske.
  • Tilpasset observasjoner i nedb?rsfelt i Birkenes kommune, Aust-Agder.
  • Flere publikasjoner p? 70-, 80- og 90-tallet.
  • Kjemisk modell koblet til hydrologisk modell.

Hydrologisk delmodell: Symboler og notasjon

Hydrologisk delmodell, detaljert introduksjon:

  • Nedb?rsfelt modellert som reservoarer $\mathbf{A}$ og $\mathbf{B}$
      • Her benyttes fet font for "merkelappene" p? reservoarene.
  • For et reservoar, for eksempel merket $\mathbf{M}$:
    • Vannstand som funksjon av tid: $M(t)$
    • Str?m ut: $Q_M(t)$
    • Evaporasjon: $E_M(t)$
    • Nedb?r ("precipitation") $P_M(t)$
    • Generelt: subindeks indikerer reservoar.
  • Avrenning til bekk. Fysiske ma?linger gjort her.

Det kan l?nne seg ? skrive en liste over st?rrelsene som inng?r i problemet.

Reservoar $\mathbf{A}$

Vi beskriver f?rst modellen for reservoar A:

Vannstand som funksjon av tid: $A(t)$ m?les i $mm$.

Vann str?mmer inn i reservoaret kun gjennom nedb?ren $P(t)$.

Vann str?mmer ut fra reservoaret gjennom evaporasjon $E_A$ og en bekk $Q_A$.

Vannstr?m/fluks m?les i $mm/dag = mm \cdot dag^{-1}$.

Reservoar $\mathbf{A}$

Regler for reservoar $\mathbf{A}$:
  • Ingen avrenning ($Q_A = 0$) n?r vannstanden $A(t)$ er mindre enn $A_{min}$
  • Avrenning intensiveres n?r vannstanden $A(t)$ ?ker over $A_{min}$
  • Evaporasjon er en funksjon av temperatur $E_A := E_A(T)$

Vi skal uttrykke dette matematisk for ? utforme modellen.

Reservoar $\mathbf{A}$

Vannstanden i $\mathbf{A}$ p?virkes positivt (?kes) av nedb?ren $P$, og negativt av evaporasjon $E$ og utstr?mming til bekken $Q_A$, alts? er

$A := A(t, P_A, E_A, Q_A)$.

Vi betrakter n? en endring i vannstanden over et lite tidsintervall $\delta t$ i magasin $\mathbf{A}$, angitt ved $\delta A$:

$\delta A = P \delta t - E_A \delta t - Q_A \delta t$

Legg merke til:

  • Symbolet $\delta$ er en liten gresk delta ($\Delta$), og brukes konvensjonelt for ? beskrive endringer og differanser. Andre steder kan vi bruke $\Delta$ eller noen ganger $h$.
  • Dimensjonen til $\delta t$ er tid, vannstr?m/fluks har dimensjon masse per tid ($mm/d?gn$) og $\delta A$ har dimensjonen masse. (Mer spesifikt mener vi her masse per kvadratmeter.)
  • Sp?rsm?l til diskusjon: g?r dimensjonene i uttrykket for endringen opp? Diskuter.

Reservoar $\mathbf{A}$

Stadig gjelder
$\delta A = P \delta t - E_A \delta t - Q_A \delta t$

Anta n? at vi p? ved tiden $t_n$ kjenner vannstanden i $\mathbf{A}$, slik at $A(t_n) = 20 mm$.

Anta ogs? at vi kjenner nedb?ren $P_A(t_n)$, evaporasjonen $E_A(t_n)$ og str?mmen til bekken $Q_A(t_n)$.

Hva vil skje med vannstanden $A$ i l?pet av det kommende lille tidsintervallet $\delta t$?

$A(t_n + \delta t) = A(t_n) + \delta A(t_n) = A(t_n) + P(t_n) \delta t - E_A(t_n) \delta t - Q_A(t_n) \delta t$

En diffligning for reservoar $\mathbf{A}$

Igjen, med utgangspunkt i

$\delta A = P \delta t - E_A \delta t - Q_A \delta t$

Del begge sider p? $\delta t$:

$\frac{\delta A}{\delta t} = P(t) - E_A(A,t) - Q_A(A,t) $

St?rrelsen $\frac{\delta A}{\delta t}$ kan leses som "en liten endring i $A$ som f?lge av en liten endring i $t$". Om vi g?r fra "en liten endring" til en "infinitesimal endring" f?r vi en differensialligning:

$\frac{d}{dt}A(t) = P(t) - E_A(A, t) - Q_A(A, t) $

En numerisk l?sning av diffligningen for reservoar $\mathbf{A}$

Vi begynner som vanlig med ? diskretisere variablene, i dette tilfellet tiden

$$t \rightarrow t_i = i \cdot \Delta t$$

Vi benytter deretter en numerisk approksimasjon p? den tidsderiverte (fremoverdifferansen):

$$\frac{d}{dt}f(t) \approx D^+f(t) = \frac{f(t + \Delta t) - f(t)}{\Delta t}$$

Dette gir

$\frac{d}{dt}A(t_n) \approx D^+A(t_n) = \frac{A(t_{n+1})-A(t_n)}{\Delta t} = P(t_n) - E_A(t_n) - Q_A(t_n) $

En numerisk l?sning av diffligningen for reservoar $\mathbf{A}$

Til slutt l?ser vi den ukjente $A(t_{n+1})$ algebraisk:

$$A(t_{n+1}) = A(t_n) + P(t_n)\Delta t - E_A(t_n)\Delta t - Q_A(t_n)\Delta t $$

Eller for ? ikke skrive for mye

$$A_{n+1} = A_n + P_n\Delta t - E_{A,n}\Delta t - Q_{A,n}\Delta t $$

Dette kalles et oppdateringsskjema, som vi ogs? husker fra forrige gang.

$\rightarrow$ Live kodeeksempel

Reservoar $\mathbf{B}$

Vi beskriver n? modellen for reservoar $\mathbf{B}$:

Vannstand som funksjon av tid: $B(t)$ m?les i $mm$.

Vann str?mmer inn i reservoaret kun gjennom bekken $Q_A(t)$.

Vann str?mmer ut fra reservoaret gjennom evaporasjon $E_B$ og en bekk $Q_B$.

Reservoar $\mathbf{B}$

Regler for reservoar $\mathbf{B}$:
  • Innstr?mming fra $Q_A$ avtar n?r vannstanden i $\mathbf{B}$ tiltar. ($A_{sig}$ er en funksjon av $B(t)$ )
  • Ingen utstr?mming til bekk n?r vannstanden er under $B_{min}$.
  • $\mathbf{B}$ har i tillegg en maksgrense $B_{max}$, hvor alt overskytende vann g?r direkte i bekken $Q_{over}$
  • Evaporasjon er en funksjon av temperatur $E_B := E_B(T)$

Innstr?mming til reservoar $\mathbf{B}$

Parameteren $A_\text{sig}$ er en funksjon av $B(t)$, $B_{min}$ og $B_{max}$:

\begin{equation} \label{eq:Asig} A_\text{sig} = \left\{ \begin{array}{lcl} 1 & ~~ & B \leq B_\text{min} \\ 1- 0.25\cdot\frac{B-B_\text{min}}{B_\text{max}-B_\text{min}} & ~~ & B_\text{min} < B \leq B_\text{max} \\ 0 & ~~ & B_\text{max} < B \end{array} \right., \qquad B_\text{max}=80\;mm. \end{equation}

Er dette en form for interpolasjon?

Reservoar $\mathbf{B}$

Vannstanden i $\mathbf{B}$ p?virkes positivt (?kes) av innstr?mming fra $\mathbf{A}$ (gjennom bekken $Q_A$, og negativt av evaporasjon $E_B$ og utstr?mming til bekken $Q_B$. Alts? er

$B := B(t, P_A, E_A, Q_A, Q_{over})$.

Vi kan uttrykke en endring i vannstanden i $\mathbf{B}$ i l?pet av et lite tidsintervall $\delta t$

$\delta B = A_{sig}Q_A \delta t - E_B \delta t - Q_B \delta t - Q_{over} \delta t$

Og finner igjen en differensialligning

$\frac{d}{dt} B(t) = A_{sig}(B) Q_A - E_B - Q_B(B) - Q_{over}(B) $

Hva skjer her?

Vi ser p? de to differensialligningene vi har funnet for modellen s? langt:

$\frac{d}{dt}A(t) = P(t) - E_A(A) - Q_A(A) $

$\frac{d}{dt} B(t) = A_{sig}(B) Q_A(A) - E_B(B) - Q_B(B) - Q_{over}(B) $

Koblede ligninger

Vi ser p? de to differensialligningene vi har funnet for modellen s? langt:

$\frac{d}{dt}A(t) = P(t) - E_A(A) - \color{blue}{Q_A(A)} $

$\frac{d}{dt} B(t) = A_{sig}(B) \color{blue}{Q_A(A)} - E_B(B) - Q_B(B) - Q_{over}(B) $

Legg merke til at vi ikke kan l?se for $B$ uten ? kjenne til $A$. Dette kalles koblede differensialligninger, eller systemer av ligninger.

Numerisk l?sning av ligningssettet

Vi kan l?se numerisk for $B(t)$ p? samme m?e som for $A(t)$, og finner at ligningssettet v?rt er

\begin{align} A_{n+1} &= A_n + \cdot P_n \Delta t - \cdot E_{A,n}\Delta t - \cdot Q_{A,n} \Delta t. \\ B_{n+1} &= B_n + \cdot A_{\text{sig},n} \cdot Q_{A,n} \Delta t - \cdot E_{B,n}\Delta t - \cdot Q_{B,n} \Delta t - \cdot Q_{\text{over},n} \Delta t . \end{align}

Hvilke st?rrelser kjenner vi ved f?rste iterasjon?

Vi kan l?se numerisk for $B(t)$ p? samme m?e som for $A(t)$, og finner at ligningssettet v?rt er

\begin{align} A_{1} &= A_0 + \cdot P_0 \Delta t - \cdot E_{A,0}\Delta t - \cdot Q_{A,0} \Delta t. \\ B_{1} &= B_0 + \cdot A_{\text{sig},0} \cdot Q_{A,0} \Delta t - \cdot E_{B,0}\Delta t - \cdot Q_{B,0} \Delta t - \cdot Q_{\text{over},0} \Delta t . \end{align}

Oppsummering

  • Oversikt over pensum
  • Introduksjon av deloppgave 2
  • Detaljert gjennomgang av hydrologisk delmodell
  • Ny teori:
    • Koblede differensialligninger
    • Interpolasjon

Takk for i dag