\input{../preambles/preamble}

\title{10 -  SYSTEMER AV DIFFERENSIALLIKNINGER}

\begin{document}
\pagenumbering{gobble}

\maketitle

%husk pendel- og pendelklokke, io og lyshastigheten, samt oppdagelsen av neptun og mallory og irvine og vulcan
% husk å ende opp med matriseeksponensial som springbrett videre

I ingeniørmatematikken finnes det ikke noe 
naturlig skarpt skille mellom lineæralgebra og differensiallikninger. 
Fysiske modeller er stort sett umulige å forstå uten noen av dem. 
Denne uken skal vi smelte dem litt ordentligere sammen. 

La oss først herje litt med likningen
\[
\ddot{x}+b\dot{x}+cx=0.
\]
Vi lager noen nye variable:
\begin{align*}
z_1(t)=x(t) \\
z_2(t)=\dot{x}(t)
\end{align*}
Hvis du nå aksepterer at vi kan skrive
\[
\mathbf z(t)=
\begin{bmatrix}
z_1(t) \\
z_2(t)
\end{bmatrix}
=
\begin{bmatrix}
x(t) \\
\dot x(t)
\end{bmatrix}
\quad
\text{og}
\quad
\dot \mathbf z(t)=
\begin{bmatrix}
\dot{z}_1(t) \\
\dot{z}_2(t)
\end{bmatrix}
=
\begin{bmatrix}
\dot x(t) \\
\ddot x(t)
\end{bmatrix}
\]
går det an å skrive $\ddot{x}+b\dot{x}+cx=0$ på formen 
\[
\dot \mathbf z = A\mathbf z.
\]



\begin{oppgave}{1}
Hva blir $A$?
\end{oppgave}

\begin{oppgave}{2}
Ser du noen løsninger?\\
(Hint: Det finnes to ekvivalente teknikker. Den ene er begynt på i forrige ukes pensum, og den andre kommer du kanskje på om du husker $\dot x = ax$ og går amok og tenker ordentlig kreativt.)
\end{oppgave}

\begin{center}
\includegraphics[scale=.27]{../figurer/L1019255}
\end{center}


Likningssystemet 
\[
\dot \mathbf z = A\mathbf z
\]
kalles gjerne et \textbf{lineært og autonomt differensiallikningssystem med konstante koeffisienter},
og matrisen $A$ kan i prinsippet være hva som helst, 
ikke bare 
\[
A=
\begin{bmatrix}
0 & 1  \\
-c & -b
\end{bmatrix}
\]
slik som i oppgave 1.
De to ovennevnte løsningsteknikkene funker på 
de aller fleste systemer (med noen unntakt, som vi skal komme tilbake til),  
og den første går som følger.
La $\lambda$ være en egenverdi til $A$, 
la $\mathbf v$ være den korresponderende egenvektoren,
og la
\[
\mathbf x(t)=e^{\lambda t}\mathbf v
\]

\begin{oppgave}{3}
Regn ut $\dot \mathbf x$
og  $A \mathbf x$.
Hva ser du?
\end{oppgave}

Hvis du skjønte oppgaven over, 
vet du i prinsippet hvordan man ter seg, 
nesten uansett hva $A$ måtte være.

\begin{oppgave}{4}
Er $\mathbf x(t)=e^{\lambda t}\mathbf v$ en løsning kun for $2\times 2$-systemer?
Hvor mange løsninger burde et $n\times n$-system ha dersom det er rettferdighet i verden?
\end{oppgave}

Finn alle løsninger til 

\begin{minipage}{55mm}
\begin{oppgave}{5} 
\begin{align*}
\dot z_1=& z_1 + 2z_2 + 2z_3\\
\dot z_2=& 2z_1 + 6z_2 + 2z_3 \\
\dot z_3=& 2z_1 + 2z_2 + 6z_3
\end{align*}
\end{oppgave}
\end{minipage}
\begin{minipage}{55mm}
\begin{oppgave}{6} 
\begin{align*}
\dot z_1=& -z_2\\
\dot z_2=& z_1
\end{align*}
\end{oppgave}
\end{minipage}
\begin{minipage}{55mm}
\begin{oppgave}{7} 
\begin{align*}
\dot z_1=& z_1+z_2\\
\dot z_2=& z_2
\end{align*}
\end{oppgave}
\end{minipage}

\begin{center}
\includegraphics[scale=.27]{../figurer/L1019390}
\end{center}

Så er det dette med initialverdier.
I alle løsningene du fant over (hvis du gjorde det riktig), 
satt det noen ubestemte konstanter. 
Disse kan bestemmes ved å spesifiser et initialkrav
\[
\mathbf z(0)=\mathbf z_0
\]
der $\mathbf z_0$ er en konstant vektor.

\begin{oppgave}{8}
Regn videre på oppgave 5, 6 og 7 med initialkravene
\[
\mathbf z(0)=
\begin{bmatrix}
5\\
0\\
2
\end{bmatrix}
\quad 
\text{og}
\quad
\begin{bmatrix}
1\\
0
\end{bmatrix}
\quad
\text{og}
\quad 
\begin{bmatrix}
0 \\
1
\end{bmatrix}
\]
henholdsvis.
\end{oppgave}

Initialverdiproblemet
\begin{align*}
\dot z_1=& z_1+z_2 \hspace{12mm} z_1(0)=0\\
\dot z_2=& z_2 \hspace{20mm} z_2(0)=1
\end{align*}
 illustrerer en viktig ting. 
Det karakteristiske polynomet til en $n \times n$-matrise kan alltid spaltes i $n$ lineære faktorer,
men det finnes ikke alltid $n$ lineært uavhengige egenvektorer.
I dette tilfellet sier vi at matrisen er \textbf{defekt}.
Differensiallikningessystemet i oppgaven har flere løsninger enn de du fant, 
og derfor klarer du ikke tilfredsstille initialkravene.

I neste semester skal vi se på hva vi gjør i slike tilfeller, 
men akkurat nå lar vi det ligge,
og ser heller på en annen måte å skrive opp 
(de samme) løsningene på.

\begin{oppgave}{9}
Niglan litt på løsningen til $\dot x = ax$ og foreslå en løsning til $\dot \mathbf x = A\mathbf x$.
\end{oppgave}

\begin{center}
\includegraphics[scale=.27]{../figurer/L1019394}
\end{center}

Du gjettet helt riktig:
\[
\mathbf x=\mathbf x_0 e^{At}
\]
Men hva i all verden skal vi gjøre med $e^{At}$? 
Hva betyr det i det hele tatt? 
Det skal vi bruke resten av semesteret på å finne ut. 
Nå skal vi over til noe helt annet, 
nemlig hvordan vi ter oss når det ikke er mulig å finne den analytiske løsningen, 
som stort sett alltid er tilfellet.

Snøskoharen (Lepus americanus) er en liten nordamerikansk hareart 
som er kjent for sine voldsomme bestandssvingninger.
Bestandssvingningene forekommer med en periode på omtrent ti år, 
og skyldes en kombinasjon av faktorer,
slik som at den formerer seg fort og at den kanadiske gaupen (Lynx canadensis) 
nesten ikke spiser noe annet enn snøskohare dersom den har tilgang på dem:\\
\url{https://www.enr.gov.nt.ca/en/services/lynx/lynx-snowshoe-hare-cycle}\\
Dette var et av de første biologiske systemene som ble forsøkt modellert av matematikere:\\
\url{https://en.wikipedia.org/wiki/Lotka–Volterra_equations}\\
Utledningen av likningene går omtrent som følger.
La $x$ være harebestanden og $y$ gaupebestanden. 
Ved fravær av gaupe vil harene formere seg eksponensielt,
slik at 
\[
\dot x=ax
\]
der $a<0$. 
Dersom det finnes gaupe, 
spiser disse harer og dette fører til en nedgang i harebestanden.
Hvis vi antar at en gaupe spiser en hare omtrent hver gang et invidid fra hver art møtes, 
burde nedgangen være proporsjonal med produktet av bestandene:
\[
\dot x=ax-bxy 
\]
Likeledes dør gaupene ut i en rate som er proporsjonal med bestanden dersom det er tomt for hare,
\[
\dot y=-cy
\]
og dersom det er hare, 
fører dette til økning i gaupebestanden proporsjonal med produktet av bestandene:
\[
\dot y=-cy+dxy
\]
Alt i alt får vi systemet
\begin{align*}
\dot{x}&=ax-bxy  \\
\dot{y}&=-cy+dxy
\end{align*}
der $a>0$, $b>0$, $c>0$ og $d>0$.
Dette systemet kan ikke løses analytisk, 
men det er ikke så vanskelig å kjøre numerisk løsning.
Dersom vi har et system av differensiallikninger med initialkrav
\begin{align*}
\dot{x}&=f(x,y) \qquad x(0)=x_0 \\
\dot{y}&=g(x,y) \qquad y(0)=y_0
\end{align*}
blir Eulers eksplisitte metode
\begin{align*}
x_{n+1}&=x_n+hf(x_n,y_n) \\
y_{n+1}&=y_n+hg(x_n,y_n)
\end{align*}

\begin{minipage}{4cm}
\begin{oppgave}{10}
Prøv. 
\end{oppgave}
\end{minipage}
\begin{minipage}{10cm}
\includegraphics[scale=.11]{../figurer/M}
\end{minipage}

Her en kode som løser forrige oppgave og plotter. 
Hvis foreldrene dine er programmerere er det 
god sjanse for at du startet å programmere i python for mange år siden,
og tok oppgaven på strak arm.
Dersom foreldrene dine er gitarister er det kanskje mer sannsynlig at 
du vet hvor skruen for å justere truss roden på gitaren din er. 
Derfor skal dere få hele koden.
(Ikke kopier over i en tekstfil. 
Se på koden og prøv å forstå strukturen,
og så prøver du å skrive koden selv uten å se på min.)

\begin{lstlisting}
import numpy as np
import matplotlib.pyplot as plt

T=20
N=200
h=T/N

t=np.zeros(N+1)
x=np.zeros(N+1)
y=np.zeros(N+1)

x[0]=2
y[0]=2

for i in range(N):
    x[i+1]=x[i]+h*x[i]*(1-y[i])
    y[i+1]=y[i]+h*y[i]*(x[i]-1)
    t[i+1]=t[i]+h


plt.plot(x,y)
plt.savefig('lotkavolterra.png')
\end{lstlisting}

Man kan plotte $x$ mot $t$ og $y$ mot $t$ i samme plot, 
eller $x$ mot $y$. 
Sistnevnte kalles \textbf{faseplott}.

\begin{oppgave}{11}
Kod opp og lag tre plott:  $t$ mot $x$, $t$ mot $y$ og $x$ mot $y$.
\end{oppgave}

\begin{center}
\includegraphics[scale=.21]{../figurer/L1009394-Pano}
\end{center}

En av de mange fordelene med å kjenne til systemer av førsteordensdifferensiallikninger, 
er at det er ryddig. 
Man trenger for eksempel ikke lage egne numeriske metoder for andre ordens differensiallikninger,
man kan bare skrive om til system og bruke de metodene man har for systemer.
La oss ta et eksempel fra fysikkens verden.

En pendel som henger fra taket og dingler uten friksjon eller luftmotstand,
tilfredsstiller
\[
\ddot{\theta}+\frac{g}{l}\sin \theta=0,
\]
der $\theta$ er vinkelutslaget, 
$l$ er pendelens lengde og $g$ er tyngdeakselerasjonen. 

\begin{oppgave}{12}
Tegn opp og utled likningen fra Newtons andre lov.
(Du trenger i tillegg massen $m$,
men denne kansellerer til slutt.)
\end{oppgave}

Pendedllikningen er umulig å løse analytisk, 
men det er bare å skrive om til system og bruke  Eulers metode.

\begin{oppgave}{13}
Sett i gang. 
\end{oppgave}

En annen klassiker fra elektroteknikken er van del Pols likning:
\[
\ddot{x}+\mu(x^2-1)\dot x + x=0
\]
Denne ble oppdaget av Balthasar van der Pol i forbindelse med radiorør, 
men modellerer også fint sprekker mellom kontinentalplater og nevronavfyring. 

\begin{oppgave}{14}
Løs van der Pol numerisk med $\mu=2$, 
og plot $x$ mot $\dot{x}$.
\end{oppgave}


\begin{center}
\includegraphics[scale=.27]{../figurer/L1019533}
\end{center}


\section*{UKENS REGULERING}


Et stivt legeme som roterer, 
tilfredsstiller eulerlikningene
\begin{align*}
I_1\dot\omega_1=(I_3-I_2)\omega_2\omega_3\\
I_2\dot\omega_2=(I_1-I_3)\omega_1\omega_3\\
I_3\dot\omega_3=(I_2-I_1)\omega_2\omega_1\\
\end{align*}
der $\omega_k$
er vinkelhastigheten rundt $k$-aksen i koordinatsystemet,
og $I_k$ er rotasjonstregheten rundt denne aksen:\\
\url{https://en.wikipedia.org/wiki/Tennis_racket_theorem}

Løs numerisk med $I_1 < I_2 < I_3$.


\section*{UKENS KJEMI}

Se oppgave 2.5 her: \\
\url{https://folk.ntnu.no/vegargje/Mattepilot_MTKJ/martaoppgaver.pdf}

\section*{UKENS ELEKTROTEKNIKK}



Av og til er det enkelt å sette systemer av differensiallikninger for kretser.
Finn differensiallikninger og alle mulige strømmer i kretsene under. 

\begin{oppgave}{1}
\begin{minipage}{5cm}

	\begin{circuitikz}
		\draw  (0,3) to (0,1);
		\draw  (0,3) to [L=$L_1$]  (4,3) to [R=$R$] (4,1) to  (0,1); 
		\draw  (4,3) to [L=$L_2$]  (8,3) to (8,1) to  (4,1); 
	\end{circuitikz}
\end{minipage}
\end{oppgave}

\begin{oppgave}{2}
\begin{minipage}{5cm}
\begin{center}
	\begin{circuitikz}
		\draw  (0,3) to (0,1);
		\draw  (0,3) to [L=$L_1$]  (4,3) to [R=$R_1$] (4,1) to  (0,1); 
		\draw  (4,3) to [L=$L_2$]  (8,3) to [R=$R_2$]  (8,1) to  (4,1); 
		\draw  (8,3) to [L=$L_3$]  (12,3) to (12,1) to  (8,1); 
	\end{circuitikz}
\end{center}
\end{minipage}
\end{oppgave}

\begin{oppgave}{3}
\begin{minipage}{5cm}
\begin{center}
	\begin{circuitikz}
		\draw  (0,3) to (0,1);
		\draw  (0,3) to [C=$C_1$]  (4,3) to [R=$R$] (4,1) to  (0,1); 
		\draw  (4,3) to [C=$C_2$]  (8,3) to (8,1) to  (4,1); 
	\end{circuitikz}
\end{center}
\end{minipage}
\end{oppgave}

\begin{oppgave}{4}
\begin{minipage}{5cm}
\begin{center}
	\begin{circuitikz}
		\draw  (0,3) to (0,1);
		\draw  (0,3) to [C=$C$]  (4,3) to [R=$R$] (4,1) to  (0,1); 
		\draw  (4,3) to [L=$L$]  (8,3) to (8,1) to  (4,1); 
	\end{circuitikz}
\end{center}
\end{minipage}
\end{oppgave}


\section*{UKENS NØTTER}

Lotka-Volterra-systemet har noe som kalles et førsteintegral:
\[
cx-b\log x+by-a\log y=C
\]
Dette kan utledes ved å dele likningene i systemet på hverandre og integrere. 

\begin{oppgave}{1}
Prøv.
\end{oppgave}

Symplektisk Euler ser slik ut:
\begin{align*}
x_{n+1}&=x_n+hf(x_n,y_n) \\
y_{n+1}&=y_n+hg(x_{n+1},y_n)
\end{align*}
Dette ser tilforlatelig ut, 
men er faktisk veldig artig. 

\begin{oppgave}{2}
Løs Lotka-Volterra med symplektisk Euler og sammenlikne med oppgave 2. 
En av metodene gir kvalitativt riktig resultat, 
og en av dem gir ikke. 
Ser du hvilken?
(Hint: prøv å finne ut hvilken som tilfredsstiller førsteintegralet best. Det kan være du må kjøre ganske langt i tid før du ser noen forskjell mellom metodene.)
\end{oppgave}


Trapesmetoden er gitt ved
\begin{align*}
x_{n+1}&=x_n+h\left(\frac{f(x_n,y_n)+f(x_{n+1},y_{n+1})}{2} \right) \\
y_{n+1}&=y_n+h\left(\frac{g(x_n,y_n)+g(x_{n+1},y_{n+1})}{2} \right)
\end{align*}

\begin{oppgave}{3}
Løs
\begin{align*}
\dot z_1=& \hspace{3mm} 2z_1-z_2 \hspace{21.5mm} z_1(0)=1\\
\dot z_2=& -z_1+2z_2 \hspace{20mm} z_2(0)=1
\end{align*}
med både Eulers metode og trapesmetoden, 
og sammenlikne med den analytiske løsningen i $t=1$.
Hva ser du?
\end{oppgave}

Det går an å kjøre trapesmetoden på ikkelineære systemer, 
men da må du vite hvordan du løser ikkelineære likningssystemer 
med flere ukjente numerisk. 
En klassier er Newtons metode: \\
\url{}\\

\begin{oppgave}{4}
Prøv van der Pol, 
Lotka-Volterra og pendelen 
med trapesmetoden. 
\end{oppgave}





\end{document}
