Skip to content
Snippets Groups Projects
Commit 8242e406 authored by jakle's avatar jakle
Browse files

hard linked to Demos

parent 71c5be59
No related branches found
No related tags found
No related merge requests found
Showing
with 10234 additions and 1 deletion
/home/jakle/Python/pilot-project-for-01005/Demos
\ No newline at end of file
This diff is collapsed.
%% Cell type:markdown id:14b14fd2 tags:
# Diagonalisering ved similartransformation
Demo af Christian Mikkelstrup og Hans Henrik Hermansen
%% Cell type:code id:aad4b143 tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id:9169d043 tags:
I denne SymPy demo vil der blive gennemgået fem eksempler på mulig diagonalisering. Til sidst præsenteres en metode til undersøgelse om matricer er similære.
%% Cell type:markdown id:df2c8fb7 tags:
### Eksempel 1
Vi betragter matricen:
%% Cell type:code id:d854a3a7 tags:
``` python
A = Matrix([[1,1,0],[2,1,-1],[0,2,1]])
A
```
%% Cell type:markdown id:24c7c5c9 tags:
Som vi finder egenværdierne af, ved "simuleret håndregning",
%% Cell type:code id:f4f95cd1 tags:
``` python
def K(l):
return A - l*eye(3)
lamda = symbols('lambda')
K(lamda)
```
%% Cell type:code id:05b970e0 tags:
``` python
karakpol = K(lamda).det()
roots(karakpol)
```
%% Cell type:markdown id:9e503221 tags:
Vi kan altså se, at der kun er én egenværdi, og at den har algebraisk multiplicitet tre! Dette eftertjekkes ved
%% Cell type:code id:cf76869f tags:
``` python
factor(karakpol)
```
%% Cell type:markdown id:d4aaa240 tags:
Egenvektorene hørende til egenværdien $1$ findes som løsningerne til den homogene matrixligning $K(1)\,x=\mathbf{0}$, altså ved
%% Cell type:code id:fd404051 tags:
``` python
K(1).nullspace()
```
%% Cell type:markdown id:25056f19 tags:
Vi har altså at
\begin{equation}
E_1=\text{span}\{(1/2,0,1)\}
\end{equation}
og dermed at summen af de geometriske multipliciteter kun er 1! Da $1<n=3$ kan $A$ *ikke* diagonaliseres!
%% Cell type:markdown id:3b268d27 tags:
### Eksempel 2
Vi betragter nu matricen
%% Cell type:code id:2405412f tags:
``` python
B = Matrix([[2,1,0],[0,1,-1],[0,2,4]])
B
```
%% Cell type:markdown id:6e6cea3c tags:
Vi benytter nu en stærkere kommando til at finde egenværdierne
%% Cell type:code id:6c05b0d3 tags:
``` python
B.eigenvals()
```
%% Cell type:markdown id:b868632e tags:
Hvor vi altså kan se, at egenværdien $2$ har algebraisk multiplicitet to, og at egenværdien $3$ har algebraisk multiplicitet en.
Nu kan vi så finde egenvektorene ved
%% Cell type:code id:7dc805e4 tags:
``` python
B.eigenvects()
```
%% Cell type:markdown id:341999c9 tags:
Altså får vi at
\begin{equation}
\begin{aligned}
E_2 = & \text{span}\,\{(1,0,0)\}\\
E_3 = & \text{span}\,\{(-1/2,-1/2,1)\}.
\end{aligned}
\end{equation}
Her ser vi, at både egenværdien $2$ og egenværdien $3$ har algebraisk multiplicitet 1. Da bare en af vores egenværdier ikke har $am(\lambda)=gm(\lambda)$ kan $B$ altså *ikke* diagonaliseres!
Dette ses også tydeligt, hvis vi benytter den indbyggede funktion:
%% Cell type:code id:1a8aba5d tags:
``` python
B.is_diagonalizable()
```
%% Cell type:markdown id:4c656c5d tags:
### Eksempel 3
Vi betragter matricen
%% Cell type:code id:603eaaed tags:
``` python
C = Matrix([[0,1,0],[-2,3,0],[0,0,2]])
C
```
%% Cell type:markdown id:ea3af290 tags:
Vi finder egenværdierne og egenvektorene direkte ved
%% Cell type:code id:d1199c56 tags:
``` python
C.eigenvects()
```
%% Cell type:markdown id:796b39a1 tags:
Her ses, at egenværdien $1$ har algebraisk multiplicitet en, og en tilhørende egenvektor, samt at egenværdien $2$ har algebraisk multiplicitet på to med to tilhørende egenvektorer. Da summen af geometrisk og algebraisk multiplicitet stemmer overens, kan vi i dette tilfælde godt diagonalisere vores matrix. De egenrum der opstår kan beskrives ved
\begin{equation}
\begin{aligned}
E_1 = & \text{span}\,\{(1,1,0)\}\\
E_2 = & \text{span}\,\{(1/2,1,0),(0,0,1)\}.
\end{aligned}
\end{equation}
Disse egenrum kan visualiseres ved
![](Billeder/Demo-E-Uge10/Eksempel3.PNG)
%% Cell type:markdown id:b9c58551 tags:
Vi fortsætter med diagonaliseringen. Vi vælger
%% Cell type:code id:1a6ae722 tags:
``` python
v1 = Matrix([1,1,0])
v2 = Matrix([S(1)/2,1,0])
v3 = Matrix([0,0,1])
Matrix.hstack(v1,v2,v3)
```
%% Cell type:markdown id:d9af2c3e tags:
Det ses at $v_1$ er en basis for $E_1$, samt at $v_2$ og $v_3$ er en basis for $E_2$. Af sætning 13.11 punkt 1, må $V=(v_1,v_2,v_3)$ altså udgøre et lineært uafhængigt sæt, som kan bruges til en basis bestående af egenvektorer.
%% Cell type:code id:da8f66aa tags:
``` python
V = Matrix.hstack(v1,v2,v3)
V
```
%% Cell type:markdown id:7776386f tags:
Og af Hovedsætning 13.14 punkt 2, får vi
%% Cell type:code id:ac99e18e tags:
``` python
Lamda = Matrix.diag([1,2,2])
Lamda
```
%% Cell type:markdown id:d8de5568 tags:
Dette kan eftertjekkes med sætning 14.7, der siger at $V^{-1}\,C\,V=\Lambda$:
%% Cell type:code id:72684605 tags:
``` python
V.inv()*C*V
```
%% Cell type:markdown id:5e7cf53f tags:
Som forventet!
%% Cell type:markdown id:8b59e217 tags:
### Eksempel 4
Vi betragter matricen
%% Cell type:code id:4de3bf3b tags:
``` python
F = Matrix([[0,-1,0],[4,5,2],[0,0,2]])
F
```
%% Cell type:markdown id:53c6b23d tags:
Vi kan direkte se, at denne kan diagonaliseres ved
%% Cell type:code id:7282d192 tags:
``` python
F.is_diagonalizable()
```
%% Cell type:markdown id:f4b7df5a tags:
Hvorved vi kan benytte den indbyggede funktion til at finde diagonaliseringen,
%% Cell type:code id:72ddff90 tags:
``` python
F.diagonalize()
```
%% Cell type:markdown id:63cfffcb tags:
Vi aflæser at egenværdierne er $1$, $2$ og $4$ med algebraisk multiplicitet (og geometrisk), samt at de tilhørende egenrum kan beskrives ved
\begin{equation}
\begin{aligned}
E_1 = & \text{span}\{(-1,1,0)\}\\
E_2 = & \text{span}\{(1,-2,1)\}\\
E_4 = & \text{span}\{(-1,4,0)\}.
\end{aligned}
\end{equation}
Disse egenrum kan visualiseres ved
![](Billeder/Demo-E-Uge10/Eksempel4.PNG)
%% Cell type:markdown id:d3fa61c2 tags:
Dette eftertjekkes med sætning 14.7
%% Cell type:code id:e3c6136a tags:
``` python
V = F.diagonalize()[0]
V.inv()*F*V
```
%% Cell type:markdown id:040bed79 tags:
Som forventet!
%% Cell type:markdown id:d84216f2 tags:
### Eksempel 5
Vi betragter den sidste matrix
%% Cell type:code id:6a6ec755 tags:
``` python
G = Matrix([[2,0,-4],[0,1,1],[0,-1,1]])
G
```
%% Cell type:markdown id:37dc7767 tags:
Som er speciel ved, at egenværdierne ikke er reelle (den kan altså ikke diagonaliseres reelt), men stadig at den (i følge SymPy) godt kan diagonaliseres:
%% Cell type:code id:c29647bc tags:
``` python
display(G.eigenvals(), G.is_diagonalizable())
```
%% Cell type:markdown id:248e91cf tags:
Vi finder egenvektorene, og benytter dem til at opskrive en kompleks basis,
%% Cell type:code id:119780da tags:
``` python
G.eigenvects()
```
%% Cell type:code id:4ab7c76c tags:
``` python
v1 = Matrix([1,0,0])
v2 = Matrix([2-2*I,I,1])
v3 = Matrix([2+2*I,-I,1])
V = Matrix.hstack(v1,v2,v3)
V
```
%% Cell type:markdown id:c2661051 tags:
samt den tilhørende diagonalmatrix:
%% Cell type:code id:f5a913b0 tags:
``` python
Lamb = Matrix.diag([2,1-I,1+I])
Lamb
```
%% Cell type:markdown id:2645d6c7 tags:
Bemærk at hvis man ikke gider skrive egenvektorerne ind i hånden, som vi lige har gjort, så kan de trækkes ud at `G.eigenvects()` på følgende måde:
%% Cell type:code id:df81f74e tags:
``` python
v1, v2, v3 = tuple([vec[2][0] for vec in G.eigenvects()]) # virker kun når gm=1
V = Matrix.hstack(v1,v2,v3)
V
```
%% Cell type:markdown id:fff8e36d tags:
Vi bemærker at der var altså slet ikke noget farligt ved at det var komplekse egenværdier. Dette ses også tydeligt, da den indbyggede funktion sagtens kan håndtere dette!
%% Cell type:code id:ead5ae14 tags:
``` python
G.diagonalize()
```
%% Cell type:markdown id:3532579a tags:
Denne kan også bruges til at opskrive $V$ og $\Lambda$ direkte:
%% Cell type:code id:2c744e2f tags:
``` python
V, Lamd = G.diagonalize()
V, Lamd
```
%% Cell type:markdown id:6fea915c tags:
## Similære matricer
Til sidst følger et eksempel, hvor vi vil finde om matricerne $A$ og $B$ er similære,
%% Cell type:code id:b4c1d040 tags:
``` python
A = Matrix([[3,2,1],[1,4,1],[1,2,3]])
B = Matrix([[6,14,13],[0,2,0],[0,0,2]])
A, B
```
%% Cell type:markdown id:ec6f9ea0 tags:
Det vides fra sætning 14.5, at de to matricer ikke kan være similære, hvis de ikke har samme karakteristiske polynomium. Dette undersøges hurtigt ved
%% Cell type:code id:0b9b7df7 tags:
``` python
A.charpoly() == B.charpoly()
```
%% Cell type:markdown id:2007df58 tags:
De har altså det samme karakteristiske polynomium, så de må altså have samme egenværdier med samme tilhørende algebraiske multiplicitet!
%% Cell type:code id:9d567aa1 tags:
``` python
A.eigenvals(), B.eigenvals()
```
%% Cell type:markdown id:9c6a1ef6 tags:
Det sidste vi mangler at undersøge er, at de to matricer har samme geometrisk multiplicitet for alle egenværdier (sætning 14.5),
%% Cell type:code id:f8cc2dd5 tags:
``` python
A.eigenvects(), B.eigenvects()
```
%% Cell type:markdown id:fde2506a tags:
Her ses at den geometriske multiplicitet også stemmer overens for alle egenværdierne, hvorved vi kan konkludere at $A$ og $B$ er similære! Hvad vi egentligt har gjort er at vise at både $A$ og $B$ er similære med diagonalmatricen
%% Cell type:code id:c0622b1f tags:
``` python
Lamb = Matrix.diag([2,2,6])
Lamb
```
%% Cell type:markdown id:64e59864 tags:
Men af sætning 14.3 får vi at $A$ og $B$ dermed også er similære!
%% Cell type:markdown id:073a2c54 tags:
Da vi nu ved at de er similære, ved vi fra sætning 14.1, at der findes en regulær matrix $M$ der opfylder $B=M^{-1}\,A\,M$.
For at finde $M$ opstiller vi koordinatmatricerne $V$ og $U$ bestående af de lineært uafhængige vektorer for $A$ og $B$ henholdsvist, begge med de to første vektorer tilhørende $E_2$ og den sidste tilhørende $E_6$.
%% Cell type:code id:c72127b0 tags:
``` python
v1 = Matrix([-2,1,0])
v2 = Matrix([-1,0,1])
v3 = Matrix([1,1,1])
V = Matrix.hstack(v1,v2,v3)
V
```
%% Cell type:code id:0ab4bd66 tags:
``` python
u1 = Matrix([-S(7)/2,1,0])
u2 = Matrix([-S(13)/4,0,1])
u3 = Matrix([1,0,0])
U = Matrix.hstack(u1,u2,u3)
U
```
%% Cell type:markdown id:9ae6a9ce tags:
Af sætning 14.7 kan vi altså opskrive
\begin{equation}
\begin{aligned}
V^{-1}\,A\,V= & \Lambda\\
U^{-1}\,A\,U= & \Lambda,
\end{aligned}
\end{equation}
hvorved vi kan sætte dem lig hinanden og regne,
\begin{equation}
\begin{aligned}
V^{-1}\,A\,V= & U^{-1}\,A\,U\\
\Leftrightarrow &\\
B = & U\,V^{-1}\,A\,V\,U^{-1}\\
\Leftrightarrow &\\
B = & (V\,U^{-1})^{-1}\,A\,(V\,U^{-1}).
\end{aligned}
\end{equation}
Matricen vi leder efter opfylder altså netop
\begin{equation}
\begin{aligned}
M = & V\,U^{-1}\\
B = & M^{-1}\,A\,M
\end{aligned}
\end{equation}
%% Cell type:code id:40e0d8c8 tags:
``` python
M=V*U.inv()
M
```
%% Cell type:markdown id:ac2a4579 tags:
Som eftertjekkes ved
%% Cell type:code id:0c117da8 tags:
``` python
B, M.inv()*A*M
```
%% Cell type:markdown id:8258b920 tags:
Som forventet
%% Cell type:markdown id:fba9ebef tags:
**Fortolkning**:
Hvad betyder det overhovedet at de er similære? Hvis vi opfatter $A$ som en afbildningsmatrix for en lineær afbildning $f:\mathbb{R}^3\to \mathbb{R}^3$ med hensyn til standard basis, kan man forestille sig $M$ som en basisskiftematrix med basisvektorer $(m_1,m_2,m_3)=((1,1,1),(1.5,4.5,3.5),(2.25,3.25,4.25))$. Altså må $B$ være den samme afbildningsmatrice $_mF_m$ for $f$ med den nye $m$-basis!
%% Cell type:code id:81c2fd3b tags:
``` python
```
This diff is collapsed.
%% Cell type:markdown id:14b14fd2 tags:
# Diagonalisering af symmetriske matricer
Demo af Christian Mikkelstrup og Hans Henrik Hermansen
%% Cell type:code id:aad4b143 tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id:9339eb92 tags:
I den sidste demo diagonaliserede vi en symmetrisk matrix ved at lave alle steps ind imellem. Dog er alle steps præcis de samme, uanset hvilken symmetrisk matrix man får! Det er derfor muligt at lave en funktion (vores, og altså ikke en indbygget) som kan give både $Q$ og $\Lambda$ som opfylder,
\begin{equation}
Q^\top\, A\, Q = \Lambda.
\end{equation}
Der er tilføjet et par ekstra tjek således at den giver en god fejlbesked, hvis den bruges på matricer som den ikke burde bruges på! Dette skal lige så meget bruges som et eksempel på hvordan man sagtens kan bygge sine egne funktioner der regner matematik ved brug af Python og SymPy!
%% Cell type:code id:c56aba08 tags:
``` python
def orto_diag_symmetrisk(A, numeric=False):
"""
Inputs:
A: Matrix()
Symmetrisk matrix af størrelse (n x n)
numeric: Bool
True/False, om der skal regnes numerisk (hurtigere)
Output:
(Q, Lamda)
Q og Lamda matricer der opfylder Q.T*A*Q = Lamda
Q er positiv ortogonal og Lamda er diagonal
"""
# Finder ud af om korrekt matrix er input
assert type(A) == type(Matrix([])), "Fejl: Input var ikke en matrix!"
assert A.shape[0] == A.shape[1], "Fejl: Matricen er ikke kvadratisk!"
assert A.T == A, "Fejl: Matricen givet er ikke symmetrisk!"
assert A.rank() == A.shape[0], "Fejl: Matricen har ikke fuld rank. \
Det er derfor ikke muligt at finde nok lineært uafhængige vektorer!"
# Nu finder vi egenvektorer
U, Lamda = A.diagonalize()
u_list = [U.col(k) for k in range(A.shape[0])]
if(numeric):
# evaluate as floating (evalf)
u_list = [u.evalf() for u in u_list]
v_list = GramSchmidt(u_list, orthonormal=True)
Q = Matrix.hstack(*v_list)
# Numeriske problemer for nogle matricer. Tjekker imaginærdelen af tallet
assert im(det(Q)) < 1e-15, "Fejl: Determinanten af Q var forventet til \
enten 1 eller -1, men imaginær del blev fundet!"
# Hvis negativ ortogonal
if(re(det(Q)) < 0):
v_list[0] = -v_list[0]
Q = Matrix.hstack(*v_list)
return Q, Lamda
```
%% Cell type:markdown id:d5326907 tags:
Denne funktion burde altså virke på alle matricer! Vi kan jo lige prøve på matricen fra sidste uges demo:
%% Cell type:code id:b0341ef7 tags:
``` python
A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
A
U, Lamda = A.diagonalize()
u_list = [U.col(k) for k in range(A.shape[0])]
u_list
```
%% Cell type:code id:bfb8aa6c tags:
``` python
Q, Lambda = orto_diag_symmetrisk(A)
Q, Lambda
```
%% Cell type:markdown id:f5797e92 tags:
Som eftertjekkes ved
%% Cell type:code id:6deff32c tags:
``` python
Q.T*A*Q, Q*Q.T
```
%% Cell type:markdown id:3654cfbf tags:
Dette virker også med større matricer,
%% Cell type:code id:6a627132 tags:
``` python
B = Matrix([[0,1,1,-1],[1,0,-1,1],[1,-1,0,1],[-1,1,1,0]])
B
```
%% Cell type:code id:dfb71f62 tags:
``` python
Q, Lambda = orto_diag_symmetrisk(B)
Q, Lambda
```
%% Cell type:markdown id:49733024 tags:
Som eftertjekkes ved
%% Cell type:code id:599c55f0 tags:
``` python
Q.T*B*Q
```
%% Cell type:markdown id:1e991bc2 tags:
Denne metode skal dog ikke bruges hovedløst! Den finder ikke frem til et svar inden for tiden af min tålmodighed, hvis man giver den følgende matrix:
%% Cell type:code id:e9986b0a tags:
``` python
C = Matrix([[1,1,0],[1,0,1],[0,1,0]])
C
# orto_diag_symmetrisk(C)
```
%% Cell type:markdown id:6558f125 tags:
Dog kan man komme i mål med en numerisk løsning ved
%% Cell type:code id:88d37e19 tags:
``` python
Q, Lamda = orto_diag_symmetrisk(C, numeric=True)
Q, Lamda
```
%% Cell type:markdown id:51159e26 tags:
Bemærk at $Q$ nu kun er approksimativt ortogonal:
%% Cell type:code id:dd9edf28 tags:
``` python
N(Q.T * Q)
```
%% Cell type:markdown id:67b9ad36 tags:
Eller vi kan tænke os endnu mere om. Den symmetriske matrix har nemlig 3 forskellige egenværdier:
%% Cell type:code id:3ac4efa5 tags:
``` python
eigs = C.eigenvals(multiple=True)
display(eigs) # eksakte egenværdier
[N(ev) for ev in eigs] # decimaltal approx
```
%% Cell type:markdown id:39b1c573 tags:
De tre egenvektorer hørende til disse tre egenværdier er derfor alle allerede ortogonale, og vi behøver kun normalisere dem:
%% Cell type:code id:896db422 tags:
``` python
v1, v2, v3 = tuple([vec[2][0].normalized() for vec in C.eigenvects()]) # virker kun når gm=1
V = Matrix.hstack(-v1,v2,v3)
V
# I en samlet kommando:
# C.diagonalize(normalize=True)
```
%% Cell type:markdown id:81e50cad tags:
Kapow, SymPy kastede matematik ud på flere A4-sider! Det er defor dog nok alligevel mere nyttigt med den numeriske approksimation, som vi kan genfinde ved:
%% Cell type:code id:246abb64 tags:
``` python
N(V)
```
%% Cell type:markdown id:ed9cdbfa tags:
De fundne $Q$ og $V$ stemmer godt overens:
%% Cell type:code id:68c6ef84 tags:
``` python
N(V-Q)
```
%% Cell type:markdown id: tags:
# Illustration af lineære transformationer med symmetriske matricer i $\mathbb{R}^2$
%% Cell type:code id: tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id: tags:
## Nogle plot-funktioner (spring gerne dette over)
%% Cell type:markdown id: tags:
Først definerer vi nogle plot-kommandoer vha. Matplotlib. Disse behøver I ikke kunne forstå.
%% Cell type:code id: tags:
``` python
# Inspired by https://scipython.com/book2/chapter-6-numpy/examples/visualizing-linear-transformations/
import matplotlib.pyplot as plt
XMIN, XMAX, YMIN, YMAX = -10, 10, -10, 10
def plot_box(B, color='k'):
"""Input: a 2x2 matrix B=[b1|b2]. Plots the convex hull of the column vectors 0,b1,b2,b1+b2."""
ix, iy = B[:,0]
jx, jy = B[:,1]
plt.plot([0, ix, ix+jx, jx, 0], [0, iy, iy+jy, jy, 0], color)
plt.axis('square') # skalerer akserne lige
plt.xlim(XMIN, XMAX) # bestem x-aksen
plt.ylim(YMIN, YMAX) # bestem y-aksen
def plot_vector(v, color='k', lw=1):
"""Plot vector v as a line with a specified color and linewidth."""
plt.plot([0, v[0]], [0, v[1]], c=color, lw=lw)
```
%% Cell type:markdown id: tags:
Vi vælger nedenfor to vektorer, $\pmb{v}_1 = [2,0]^T$ og $\pmb{v}_1 = [0,3]^T$, og plotter randen af firkanten $F$ givet ved:
$$
F = \Big\{ c_1 \pmb{v}_1 + c_2 \pmb{v}_2 \; \Big\vert \; 0 \le c_1,c_2 \le 1 \Big\}
$$
%% Cell type:code id: tags:
``` python
v_1 = Matrix(2,1,[2,0])
v_2 = Matrix(2,1,[0,3])
V = v_1.row_join(v_2) # V = [v_1|v_2]
plot_vector(v_1,color='r',lw=2) # plot v_1 med rød
plot_vector(v_2,color='b',lw=2) # plot v_2 med gul
plot_box(V,color='y') # plot randen af firkanten F med gul
```
%% Cell type:markdown id: tags:
Vi betragter en linear transformation $\mathbb{R}^2 \to \mathbb{R}^2$ givet ved $\pmb{x} \mapsto A\pmb{x}$, hvor $A$ er:
%% Cell type:code id: tags:
``` python
A = Matrix(2,2,[2,1,1,2])
A
```
%% Cell type:markdown id: tags:
Vi diagonaliserer den symmetriske matrix $A$:
%% Cell type:code id: tags:
``` python
Q, Lamda = A.diagonalize(normalize=True)
Q, Lamda
```
%% Cell type:markdown id: tags:
Vi undersøger billedet $A(F)$, der er givet ved:
$$
A(F) = \Big\{ c_1 A \pmb{v}_1 + c_2 A \pmb{v}_2 \; \Big\vert \; 0 \le c_1,c_2 \le 1 \Big\}
$$
%% Cell type:code id: tags:
``` python
plot_box(V,color='k') # plot randen af F med sort
plot_box(A*V,color='y') # plot A(F) med gul
plot_vector(A*v_1,color='r') # plot A(v_1) med rød
plot_vector(A*v_2,color='b') # plot A(v_2) med blå
```
%% Cell type:markdown id: tags:
Da $A = Q \Lambda Q^T$ kan vi undersøge hvad der sker trin for trin. Først $Q^T(F)$, der roterer mod venstre i skiftet fra standard $e$-basis til $q$-basis (husk $Q = {}_eM_q$ og $Q^T = Q^{-1} = {}_qM_e$)
%% Cell type:code id: tags:
``` python
plot_box(V,color='k')
plot_box(Transpose(Q)*V,color='y')
plot_vector(Transpose(Q)*v_1,color='r')
plot_vector(Transpose(Q)*v_2,color='b')
```
%% Cell type:markdown id: tags:
Så $\Lambda Q^T(F)$. Husk at $\Lambda$ afbildningsmatricen mht til $q$-basis. Diagonalmatricen $\Lambda = \mathrm{diag}(1,3)$ strækker $y$-koordinaten med en faktor $3$, mens den lader $x$-koordinaten være uændret.
%% Cell type:code id: tags:
``` python
plot_box(V,color='k')
plot_box(Lamda*Transpose(Q)*V,color='y')
plot_vector(Lamda*Transpose(Q)*v_1,color='r')
plot_vector(Lamda*Transpose(Q)*v_2,color='b')
```
%% Cell type:markdown id: tags:
Og endelig $Q \Lambda Q^T(F)$ som jo selvfølgelig er det samme som $A(F)$. Multiplikation med $Q$ roterer mod højre i skiftet fra $q$-basis tilbage til $e$-basis.
%% Cell type:code id: tags:
``` python
plot_box(V,color='k')
plot_box(Q*Lamda*Transpose(Q)*V,color='y')
plot_vector(Q*Lamda*Transpose(Q)*v_1,color='r')
plot_vector(Q*Lamda*Transpose(Q)*v_2,color='b')
```
%% Cell type:markdown id: tags:
## Mere advancerede matplotlib-kommandoer
%% Cell type:markdown id: tags:
I matplotlib-pakken kan kan vælge følgende pre-definerede "plot styles":
%% Cell type:code id: tags:
``` python
print(plt.style.available)
```
%% Cell type:markdown id: tags:
Vi vælger en, tilføjer plottet lidt mere information og gemmer det som en png-fil:
%% Cell type:code id: tags:
``` python
plt.style.use('ggplot')
plt.title('Plot af firkanten F og dens billede A(F) under A')
plt.xlabel('x')
plt.ylabel('y')
plt.xticks(list(range(-10,12,2)))
plt.yticks(list(range(-10,12,2)))
plot_box(V,color='k')
plot_box(Q*Lamda*Transpose(Q)*V,color='y')
plot_vector(Q*Lamda*Transpose(Q)*v_1,color='r')
plot_vector(Q*Lamda*Transpose(Q)*v_2,color='b')
# Gemmer plot som png-fil:
plt.savefig('mit_plot.png')
```
%% Cell type:markdown id: tags:
Lidt mindre seriøst, har matplotlib en [xkcd](https://xkcd.com/) sketch-style drawing mode. Det fungerer bedst med Humor fonts installeret.
%% Cell type:code id: tags:
``` python
plt.xkcd()
plt.style.use('seaborn-v0_8-colorblind')
plt.title('Plot af firkanten F og dens billede A(F) under A')
plt.xlabel('x')
plt.ylabel('y')
plt.xticks(list(range(-10,12,2)))
plt.yticks(list(range(-10,12,2)))
plt.grid()
plot_box(V,color='k')
plot_box(Q*Lamda*Transpose(Q)*V,color='y')
plot_vector(Q*Lamda*Transpose(Q)*v_1,color='r')
plot_vector(Q*Lamda*Transpose(Q)*v_2,color='b')
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id:14b14fd2 tags:
# Symmetriske matricer
Demo af Christian Mikkelstrup, Hans Henrik Hermansen og Jakob Lemvig
%% Cell type:code id:aad4b143 tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id:1a37dcc0 tags:
I denne demo gennemgåes eksempler der udregner og benytter ortonormale baser, samt diagonalisering af en symmetriske matrix ved ortogonal substitution.
%% Cell type:markdown id:0e7c596c tags:
## Ortonormale baser
I de følgende tre eksempler benyttes tre vektorer $u_1,u_2,u_3\in \mathbb{R}^4$,
%% Cell type:code id:1ec4a432 tags:
``` python
u_1 = Matrix([1,1,1,1])
u_2 = Matrix([1,0,1,0])
u_3 = Matrix([1,0,0,1])
u_1, u_2, u_3
```
%% Cell type:markdown id:1fd990e7 tags:
### Eksempel 1
I dette eksempel udregner vi længden af og vinkler mellem vektorer. Først længden af de tre vektorer:
%% Cell type:code id:fe38fd6e tags:
``` python
display(u_1.norm(), u_2.norm(), u_3.norm())
```
%% Cell type:markdown id:c9e15602 tags:
Vinkler mellem vektorer blev allerede præsenteret i Demo'en fra lille dag, uge 6. Dog, da vi skal gøre det mange gange, kan vi definere vores egen funktion:
%% Cell type:code id:e554dee4 tags:
``` python
def angle(u_1, u_2):
cos_to_v = (u_1.dot(u_2))/(u_1.norm()*u_2.norm())
return acos(cos_to_v)
```
%% Cell type:code id:4d7fc403 tags:
``` python
# vinkel mellem u_1 og u_2
angle(u_1,u_2)
```
%% Cell type:code id:3f22ed99 tags:
``` python
# vinkel mellem u_2 og u_3
angle(u_2,u_3)
```
%% Cell type:code id:b19a8657 tags:
``` python
# vinkel mellem u_3 og u_1
angle(u_3,u_1)
```
%% Cell type:markdown id:6678b63d tags:
### Eksempel 2
I dette eksempel finder vi en ortonormal basis for underrummet udspændt af $u_1, u_2$ og $u_3$, altså for $U=\text{span}\{u_1, u_2, u_3\}$.
Vi kan således beskrive koordinatmatricen for vektorene med hensyn til den sædvanlige basis som,
%% Cell type:code id:cd96f608 tags:
``` python
eU = Matrix([[1,1,1],[1,0,0],[1,1,0],[1,0,1]])
eU
```
%% Cell type:markdown id:bc6f35f6 tags:
Hvor rangen af $_eU$ findes til
%% Cell type:code id:08567dff tags:
``` python
eU.rank()
```
%% Cell type:markdown id:a13706c7 tags:
Vi kan altså sige, at de tre vektorer er lineært uafhængige, og derved kan bruges som basis for $U$.
%% Cell type:markdown id:0efdfe0c tags:
En **ortonormal basis** for underrummet U findes ved brug af **Gram-Schmidt** metoden. Her vises den ved simuleret håndregning,
%% Cell type:code id:7a7502ab tags:
``` python
v_1 = u_1/u_1.norm()
v_1
```
%% Cell type:code id:b96f4a42 tags:
``` python
w_2 = u_2 - u_2.dot(v_1)*v_1
v_2 = w_2/w_2.norm()
v_2
```
%% Cell type:code id:b7cd2536 tags:
``` python
w_3 = u_3 - u_3.dot(v_1)*v_1 - u_3.dot(v_2)*v_2
v_3 = w_3/w_3.norm()
v_3
```
%% Cell type:markdown id:69e9c289 tags:
Hvorved vi har vores basis bestående af $v_1$, $v_2$ og $v_3$.
Denne basis kan også fåes direkte med den indbyggede funktion,
%% Cell type:code id:ff0308ff tags:
``` python
GramSchmidt([u_1, u_2, u_3], orthonormal=True)
```
%% Cell type:markdown id:aee7819a tags:
### Eksempel 3
Her følger et eksempel på en ortonormal basis i $\mathbb{R}^3$. Vi får givet to vektorer, som er ortogonale og som har længden $1$,
%% Cell type:code id:9efc8a8f tags:
``` python
v_1 = Matrix([1,2,-2])/S(3)
v_2 = Matrix([2,1,2])/S(3)
v_1,v_2
```
%% Cell type:markdown id:fb6b41e5 tags:
Den tredje vektor findes ved krydsproduktet,
%% Cell type:code id:dc6d7b21 tags:
``` python
v_3 = v_1.cross(v_2)
```
%% Cell type:markdown id:a543901b tags:
Nu må det altså gælde, at $(v_1, v_2, v_3)$ udgør en en ortonormal basis for $\mathbb{R}^3$.
Koordinatmatricen $V$ for de tre vektorer er et eksempel på en ortogonal matrix,
%% Cell type:code id:9cb9d3a5 tags:
``` python
V = Matrix.hstack(v_1, v_2, v_3)
V
```
%% Cell type:markdown id:0899270f tags:
Dette kan eftertjekkes ved
%% Cell type:code id:23601557 tags:
``` python
V.T * V
```
%% Cell type:markdown id:bc93e74a tags:
Når en ortonormal basis i $\mathbb{R}^3$ er fremskaffet ved krydsproduktet, er koordinatmatricen positiv ortogonal. Altså har den valgte basis den **sædvanlige orientering**. Dette eftertjekkes ved
%% Cell type:code id:db5e6485 tags:
``` python
det(V)
```
%% Cell type:markdown id:dd14ae58 tags:
som forventet.
%% Cell type:markdown id:8a5a6ffb tags:
# Diagonalisering af en symmetrisk matrix ved ortogonal substitution
Den symmetriske matrix $A$ er givet ved
%% Cell type:code id:c43ae7ec tags:
``` python
A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
A
```
%% Cell type:markdown id:3235b1b8 tags:
Vi vil bestemme en positiv ortogonal matrix $Q$ og en diagonalmatrix $\Lambda$ så $Q^\top\, A\, Q = \Lambda$.
Vi starter med at finde 3 lineært uafhængige egenvektorer vha. SymPy:
%% Cell type:code id:65870bcf tags:
``` python
A.eigenvects()
```
%% Cell type:code id:31db740e tags:
``` python
u_1 = Matrix([-1,Rational(1,2),1])
u_2 = Matrix([Rational(1,2),1,0])
u_3 = Matrix([1,0,1])
```
%% Cell type:markdown id:0014f5ec tags:
Eller i et "go" uden at skulle skrive $u$-vektorerne ind manuelt:
%% Cell type:code id:789ade26 tags:
``` python
U, Lamda = A.diagonalize()
[u_1,u_2,u_3] = [U.col(k) for k in range(3)] # hver søjle i U tilgås og gemmes som vektor
[u_1,u_2,u_3]
```
%% Cell type:markdown id:4231d98d tags:
Nu kan vi bruge den indbyggede metode til at (eller håndregning som vist højere oppe) til at finde den ortonormale basis. Vi finder:
%% Cell type:code id:b316d2aa tags:
``` python
v_1 = u_1.normalized()
v_2, v_3 = GramSchmidt([u_2,u_3], orthonormal=True)
Q = Matrix.hstack(v_1,v_2,v_3)
Q
```
%% Cell type:markdown id:f9ae3c78 tags:
Bemærk at vi kører Gram-Schmidt på hhv. den ene vektor $u_1$ og bagefter på de to vektorer $u_2, u_3$. Kommandoen `v_1,v_2, v_3 = GramSchmidt([u_1,u_2,u_3], orthonormal=True)` giver det samme (hvorfor?).
Vi kan undersøge om $Q$ er positiv ortogonal ved
%% Cell type:code id:1a029d44 tags:
``` python
det(Q)
```
%% Cell type:markdown id:d17febf7 tags:
Det er den ikke. Dog kan dette fixes ved at skifte fortegn på én af de tre $v$-vektorer,
%% Cell type:code id:9c721543 tags:
``` python
Q = Matrix.hstack(v_1, v_2, -v_3)
Q
```
%% Cell type:markdown id:011c09f8 tags:
Vi tjekker at den opdaterede $Q$-matrix er positiv ortogonal:
%% Cell type:code id:e134d725 tags:
``` python
det(Q)
```
%% Cell type:markdown id:8ae5f018 tags:
Rækkefølgen af egenværdierne i diagonalmatricen $\Lambda$ er bestemt af rækkefølgen af egenvektorer, som vi valgte. Den er lig $U^T A U$ og givet ved:
%% Cell type:code id:a80ad5e1 tags:
``` python
Lamda
```
%% Cell type:markdown id:d2f01d15 tags:
Lad os tjekke at Gram-Schmidt ikke har ændret ved $\Lambda$. Vi skriver:
%% Cell type:code id:524e7c71 tags:
``` python
Lamda == Q.T*A*Q
```
%% Cell type:markdown id:695ebb02 tags:
Som forventet er $\Lambda = U^T A U = Q^T A Q$. Bemærk at dette "Gram-Schmidt"-trick ikke skal bruges for matricer der opfyder $AA^T \neq A^TA$, da der altid vil gælde $\Lambda = U^T A U \neq Q^T A Q$.
%% Cell type:markdown id: tags:
# Differentialligningssystemer
Demo af Christian Mikkelstrup og Hans Henrik Hermansen
%% Cell type:code id: tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id: tags:
## Homogene systemer af første ordens lineære diff-ligninger
## 1. Systemmatricen har to reelle egenværdier
Givet differentialligningsystemet
\begin{gather*}
\frac{\text{d}}{\text{d}t}x_1(t) = 5x_1(t) - 3x_2(t)\\
\frac{\text{d}}{\text{d}t}x_2(t) = 6x_1(t) - 4x_2(t)
\end{gather*}
%% Cell type:markdown id: tags:
### a. Simuleret håndregning via diagonaliseringsmetoden
Først løser vi systemet i hånden vha. diagonaliseringsmetoden, dvs. vi opstiller og analyserer systemmatricen.
%% Cell type:code id: tags:
``` python
t,C1,C2 = symbols("t C1:3")
A = Matrix([[5,-3],[6,-4]])
ev = A.eigenvects()
res = C1 * E ** (ev[0][0]*t) * ev[0][2][0] + C2 * E ** (ev[1][0]*t) * ev[1][2][0]
```
%% Cell type:markdown id: tags:
Løsningen $[x_1(t), x_2(t)]^T$ er altså givet ved:
%% Cell type:code id: tags:
``` python
res
```
%% Cell type:markdown id: tags:
hvor $C_1, C_2$ er arbitrære reelle eller komplekse konstanter$.
%% Cell type:markdown id: tags:
### b. via dsolve
Vi kan også finde løsningen direkte med dsolve
%% Cell type:code id: tags:
``` python
x1 = Function("x1")
x2 = Function("x2")
eq1 = Eq(diff(x1(t),t),5 * x1(t) - 3 * x2(t))
eq2 = Eq(diff(x2(t),t),6*x1(t) - 4*x2(t))
#
dsolve([eq1, eq2])
```
%% Cell type:markdown id: tags:
## 2. Systemmatricen har to komplekse egenværdier
Givet differentialligningssystemet
\begin{gather*}
\frac{\text{d}}{\text{d}t}x_1(t) = 2x_1(t) + 2x_2(t)\\
\frac{\text{d}}{\text{d}t}x_2(t) = -x_1(t) + 4x_2(t)
\end{gather*}
%% Cell type:markdown id: tags:
## a. Simuleret håndregning
Fremgangsmåden her er den samme.
%% Cell type:code id: tags:
``` python
A = Matrix([[2,2],[-1,4]])
ev = A.eigenvects()
res = C1 * E ** (ev[0][0]*t) * ev[0][2][0] + C2 * E ** (ev[1][0]*t) * ev[1][2][0]
res
```
%% Cell type:markdown id: tags:
hvor $C_1, C_2 \in \mathbb{C}$ hvis vi er ude efter den fuldstændige komplekse løsning.
%% Cell type:markdown id: tags:
Nu finder vi også den fuldstændige reelle løsning vha. metode 17.5
%% Cell type:code id: tags:
``` python
D1,D2 = symbols("D1:3")
a = re(ev[0][0])
b = im(ev[0][0])
re_ev = re(ev[0][2][0])
im_ev = im(ev[0][2][0])
u1 = E ** (a * t) * (cos(b * t) * re_ev - sin(b * t) * im_ev)
u2 = E ** (a * t) * (sin(b * t) * re_ev + cos(b * t) * im_ev)
D1 * u1 + D2 * u2
```
%% Cell type:markdown id: tags:
hvor $D_1, D_2 \in \mathbb{R}$. Bemærk at hvis vi i stedet skriver "hvor $D_1, D_2 \in \mathbb{C}$" fås fuldstændige komplekse løsning.
%% Cell type:markdown id: tags:
### b. Med dsovle
%% Cell type:code id: tags:
``` python
eq1 = Eq(diff(x1(t),t),2*x1(t) + 2*x2(t))
eq2 = Eq(diff(x2(t),t),-x1(t) + 4*x2(t))
dsolve([eq1, eq2])
```
%% Cell type:markdown id: tags:
Selvom løsningerne er lidt forskellige ud, så giver de anledning til samme fuldstændige løsningsmængde. For at indse dette skal vi blot sætte $D_1 = C_2$ og $D_1 = C_2$.
%% Cell type:markdown id: tags:
## 3. Betinget løsning med plot
Lad os se på systemet fra eksempel 1 igen.
\begin{gather*}
\frac{\text{d}}{\text{d}t}x_1(t) = 5x_1(t) - 3x_2(t)\\
\frac{\text{d}}{\text{d}t}x_2(t) = 6x_2(t) - 4x_2(t)
\end{gather*}
Vi kom frem til den fuldstændige løsning
\begin{gather*}
\begin{bmatrix} x_1(t) \\ x_2(t) \end{bmatrix} = C_1 e^{-t} \begin{bmatrix}1 \\ 2\end{bmatrix} + C_2 e^{2t} \begin{bmatrix} 1 \\ 1 \end{bmatrix}
\end{gather*}
Lad os finde et $C_1$ og $C_2$, der medfører at $x_1(0) = 2$ og $x_2(0) = 3$.
%% Cell type:markdown id: tags:
### a. Simuleret håndregning
Først finder vi lige løsningen igen
%% Cell type:code id: tags:
``` python
t,C1,C2 = symbols("t C1:3")
A = Matrix([[5,-3],[6,-4]])
ev = A.eigenvects()
fuld = C1 * E ** (ev[0][0]*t) * ev[0][2][0] + C2 * E ** (ev[1][0]*t) * ev[1][2][0]
x1_fuld, x2_fuld = fuld[0],fuld[1]
x1_fuld, x2_fuld
```
%% Cell type:markdown id: tags:
Vi sætter $t=0$ i løsningerne og bruger begyndelsesbetingelsen $x_1(0) = 2$ og $x_2(0) = 3$, hvilket giver ligningerne:
%% Cell type:code id: tags:
``` python
display(Eq(x1_fuld.subs(t,0),2),Eq(x2_fuld.subs(t,0),3))
```
%% Cell type:markdown id: tags:
Dette løses let, fx ved:
%% Cell type:code id: tags:
``` python
solve([Eq(x1_fuld.subs(t,0),2),Eq(x2_fuld.subs(t,0),3)])
```
%% Cell type:markdown id: tags:
Vi kunne også gøre dette direkte med $\text{dsolve}$
%% Cell type:code id: tags:
``` python
eq1 = Eq(diff(x1(t),t),5 * x1(t) - 3 * x2(t))
eq2 = Eq(diff(x2(t),t),6*x1(t) - 4*x2(t))
dsolve([eq1, eq2],ics={x1(0) : 2, x2(0) : 3})
```
%% Cell type:markdown id: tags:
Nu, hvor vi har fundet $c_1$ og $c_2$, kan vi plotte vores løsning.
%% Cell type:code id: tags:
``` python
plot(x1_fuld.subs([(C1,2),(C2,1)]),x2_fuld.subs([(C1,2),(C2,1)]),xlim=(-1,1),ylim=(0,8),legend=True)
```
%% Cell type:markdown id: tags:
# Andenordens lineære differentialligninger med konstante koefficienter
%% Cell type:code id: tags:
``` python
from sympy import *
init_printing()
t = symbols('t', real=True)
```
%% Cell type:markdown id: tags:
Figur 1: ![Fig 1](vogn_fig1.png)
Figur 2: ![Fig 2](vogn_fig2.png)
%% Cell type:markdown id: tags:
En vogn er på fig.1 fastgjort til muren med en fjeder. Fjederen er i sin hvilestilling, så vognen holder stille og roligt ved $x = 0$. Til et bestemt tidspunkt $t$ er vognen bragt ud til positionen $x(t)$ mod højre, se fig.2. Den er da i følge Hooke's lov påvirket af en kraft mod venstre af størrelsen $F = -k x(t)$, hvor $k$ er en (positiv) konstant som angiver fjederens stivhed. Jo længere mod højre vognen er placeret, desto mere trækker fjederen altså i den mod venstre. Og når vi slipper den, begynder den at køre! Vi vil gerne beskrive vognens bevægelse, det vil sige bestemme dens placering $x(t)$ til tiden $t$. Vi er nødt til at gå indirekte frem. Det viser sig at vi kan opstille et udtryk for bevægelsens acceleration - og det kan vi arbejde videre med.
%% Cell type:markdown id: tags:
# Opgave 1
## a)
I følge Newtons anden lov er den kraft der påvirker vognen, lig med vognens masse $M$ gange dens acceleration. Opstil en ligning som kæder accelerationen sammen med stedfunktionen $x(t)$ ved at udfylde `XX` i den følgende kommando:
```
x = Function("x")
k, M = symbols('k, M', positive=True)
diff_lign1 = Eq(diff(x(t),t,t), XX * x(t))
```
%% Cell type:markdown id: tags:
### Løsning a)
%% Cell type:markdown id: tags:
Vi har fra Newton og Hooke to forskellige udtryk for den resulterende kraft $F$, nemlig:
$$ F = M a(t) \quad \text{og} \quad M a(t) = -k x(t)$$
og da accelerationen er den anden afledede af stedfunktionen, får vi
%% Cell type:code id: tags:
``` python
x = Function("x")
k, M = symbols('k, M', positive=True)
diff_lign1 = Eq(diff(x(t),t,t), -(k/M) * x(t))
diff_lign1
```
%% Cell type:markdown id: tags:
Det er nu lykkedes os at opstille en andenordens differentialligning, hvori stedfunktionen optræder som ukendt funktion. Vi må forvente at vognen vil køre frem og tilbage, henover $x = 0$. Men hvad betyder egentlig fjederens stivhed $k$ for svingningens hastighed? Dette spørgsmål nærmer vi os i det følgende.
I det følgende sætter vi som eksempel vognens masse til $M=1$.
## b)
Find et foreløbigt udtryk for $x(t)$ vha. `dsolve`. Hvor mange yderligere oplysninger (begyndelsesbetingelser) skal vi bruge for endeligt at fastlægge $x(t)$, når vi forudsætter at $k$ er kendt?
## c)
Antag at en fotooptagelse viser at vognen til tiden t = 0 var ved stedet $x=1$, og til tiden $t = 1$ ved stedet $x=5$, og plot på den baggrund $x(t)$ for henholdsvis $k=1, 2, 3$. Kommentér plottet.
%% Cell type:markdown id: tags:
### Løsning b)+c)
%% Cell type:code id: tags:
``` python
fuld_sol1 = dsolve(diff_lign1.subs(M,1))
fuld_sol1
```
%% Cell type:markdown id: tags:
I den fuldstændige løsning optræder der to ukendte konstanter $C_1$ og $C_2$. Vi har derfor brug for to begyndelsesbetingelser, hvis vi forudsætter $k$ kendt.
%% Cell type:code id: tags:
``` python
part_sol1 = dsolve(diff_lign1.subs(M,1),ics = {x(0) : 1, x(1) : 5}).rhs
part_sol1
```
%% Cell type:code id: tags:
``` python
p1 = plot(part_sol1.subs(k,1),(t,0,10), show=False, xlabel='t', ylabel='x(t)', legend=True)
p1[0].label = 'k=1'
p2 = plot(part_sol1.subs(k,2),(t,0,10), show=False)
p2[0].label = 'k=2'
p3 = plot(part_sol1.subs(k,3),(t,0,10), show=False)
p3[0].label = 'k=3'
p1.append(p2[0])
p1.append(p3[0])
p1.show()
```
%% Cell type:markdown id: tags:
Vi bemærker at alle tre grafer, som forventet, opfylder begyndelsesbetingelserne idet de passerer $(0,1)$ og $(1, 5)$. Og vi ser at desto *stivere* fjederen er, desto *hurtigere* kører vognen frem og tilbage.
%% Cell type:markdown id: tags:
# Opgave 2
Nu gør vi model-situationen lidt mere realistisk. Der er selvfølgelig en vis modstand mod bevægelsen som skyldes vindmodstand, knirken i vognens hjul mv. Det er oplagt at forvente at denne gnidningsmodstand afhænger af hastigheden: desto større hastigheden er, desto større er modstanden imod bevægelsen. Vi antager enkelt at modstandskraften kan skrives som $-n$ ganget med hastigheden, hvor $n$ er et positivt tal der angiver gnidningsbetingelserne. Mon vognen stadig vil køre frem og tilbage, eller vil svingningen ophøre hvis $n$ er tilstrækkelig stor? Dette spørgsmål nærmer vi os i det følgende.
%% Cell type:markdown id: tags:
# d)
Angiv den samlede kraft F der nu virker på vognen i fig.2, og opstil den andenordens differentialligning der sammenfatter situationen.
%% Cell type:markdown id: tags:
## Løsning d)
%% Cell type:code id: tags:
``` python
n = symbols('n', positive=True)
diff_lign2 = Eq(diff(x(t),t,t), -k/M*x(t)-n/M*diff(x(t),t))
diff_lign2
```
%% Cell type:markdown id: tags:
I det følgende sætter vi som eksempel fjederens stivhed til $k=1$ (og stadig $M=1$). Lad os undersøge om vi stadig med to begyndelsesbetingelser kan opnå en endelig løsning på det nye bevægelsesproblem. En videooptagelse viser at vognen til tiden t = 0 var ved stedet $x=0$, og at dens hastighed til tiden $t = 0$ var $1$. Den fuldstændige løsning er:
%% Cell type:code id: tags:
``` python
diff_lign2 = diff_lign2.subs(M,1).subs(k,1)
fuld_sol1 = dsolve(diff_lign2)
fuld_sol1
```
%% Cell type:markdown id: tags:
Vi bemærker at $\sqrt{n-2}$ optræder i løsningen. Vi må derfor forvente at løsningsformen ændrer sig ved $n=2$.
%% Cell type:markdown id: tags:
# e)
Find den partikulære løsningen $x(t)$ for henholdvis $n=1,2,3$.
# f)
Plot de tre løsninger og kommentér plottet.
%% Cell type:markdown id: tags:
## Løsning e)+f)
%% Cell type:markdown id: tags:
Da vi i d) så at løsningen ændrer form ved $n=2$, må vi hellere løse differentialligningen for hver af de 3 $n$-værdier. For $n=1$ bliver det til:
%% Cell type:code id: tags:
``` python
dsolve(diff_lign2.subs(n,1), ics = {x(0) : 0, diff(x(t),t).subs(t,0) : 1})
```
%% Cell type:markdown id: tags:
Vi bruger en ´for´-løkke for de tre tilfælde; alternativt kan man bare bruge ovenstående `dsolve`-kommando 3 gange.
%% Cell type:code id: tags:
``` python
part_sol2 = []
for nn in [1,2,3]:
part_sol2.append(dsolve(diff_lign2.subs(n,nn), ics = {x(0) : 0, diff(x(t),t).subs(t,0) : 1}).rhs)
print('De tre løsninger:')
part_sol2
```
%% Cell type:code id: tags:
``` python
p1 = plot(part_sol2[0],(t,0,10), show=False, xlabel='t', ylabel='x(t)', legend=True)
p1[0].label = 'n=1'
p2 = plot(part_sol2[1],(t,0,10), show=False)
p2[0].label = 'n=2'
p3 = plot(part_sol2[2],(t,0,10), show=False)
p3[0].label = 'n=3'
p1.append(p2[0])
p1.append(p3[0])
p1.show()
```
%% Cell type:markdown id: tags:
Det er kun i tilfældet $n = 1$ at vognen, efter en tur ude på vejens højre side, passerer $x = 0$ og derfor udfører egentlige svingninger pga. $\sin$ i funktionsudstrykket. Men svingningen dæmpes, og i ingen af de tre tilfælde går der lang tid før at vognen (næsten) står stille ved $x=0$
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
# 2. ordens lineære differentialligning med konst. koefficenter
Demo af Christian Mikkelstrup og Hans Henrik Hermansen
%% Cell type:code id: tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id: tags:
Vi betragter følgende differentialliging
\begin{equation*}
\frac{\text{d}^2}{\text{d}t^2}x(t) - 3\frac{\text{d}}{\text{d}t}x(t) + 2x(t) = \cos(2t)
\end{equation*}
og ønsker at finde **reelle** løsninger.
%% Cell type:markdown id: tags:
## 1. Simuleret håndregning
%% Cell type:markdown id: tags:
### a. Undersøgelse af den tilsvarende homogene ligning
Først finder vi løsningen til den homogene ligning. Dette gøres ved at løse karakterligningen
%% Cell type:code id: tags:
``` python
l, C1, C2 = symbols("l C1:3")
solve(l ** 2 - 3 * l + 2)
```
%% Cell type:markdown id: tags:
Så den homogene del af løsninger er
\begin{equation*}
C_1 e^{2t} + C_2 e^{t}
\end{equation*}
%% Cell type:markdown id: tags:
hvor $C_1,C_2 \in \mathbb{R}$.
%% Cell type:markdown id: tags:
### b. Undersøgelse af den inhomogene ligning
Her bruger vi gættemetoden. Vi kan se at højresiden har formen $\cos(2t)$. Derfor gætter vi på en løsning på formen $a \sin(2t) + b \cos(2t)$
%% Cell type:code id: tags:
``` python
a,b,t = symbols("a b t")
gæt = a * sin(2*t) + b * cos(2*t)
eq_gæt = Eq(simplify(diff(gæt,t,t) - 3 * diff(gæt,t) + 2 * gæt),cos(2*t))
eq_gæt
```
%% Cell type:markdown id: tags:
Ligningen skal gælde for alle værdier af $t$. Ved at sætte $t=0$ fås $-6a-2b = 1$. Ved at sætte $t=\pi/4$ fås, da $\sin(\pi/2)=1$ og $\cos(\pi/2)=0$ fås $-2a+6b = 0$. Det ses let i SymPy:
%% Cell type:code id: tags:
``` python
eq_gæt.subs(t,0), eq_gæt.subs(t,pi/4)
```
%% Cell type:markdown id: tags:
og løsningen er:
%% Cell type:code id: tags:
``` python
ab_val = solve([eq_gæt.subs(t,0), eq_gæt.subs(t,pi/4)])
ab_val
```
%% Cell type:markdown id: tags:
Vi mangler at tjekke at disse valg af $a$ og $b$ opfylder $−2a \sin(2t)−6a\cos(2t)+6b\sin(2t)−2b\cos(2t)=cos(2t)$ for alle $t$. Det kan gøres ved:
%% Cell type:code id: tags:
``` python
eq_gæt.subs(ab_val)
```
%% Cell type:markdown id: tags:
Deraf ender vi med en fuldstændig løsning
\begin{equation*}
x(t) = C_1 e^{2t} + C_2 e^t - \frac{3}{20}\sin(2t) - \frac{1}{20}\cos(2t),
\end{equation*}
hvor $C_1,C_2 \in \mathbb{R}$.
%% Cell type:markdown id: tags:
### c. En ønsket partikulær løsning
Nu er vi interesserede i at finde en løsning med begyndelsesbetingelserne $x(0) = 0$ og $x'(0) = 1$.
Dette gør vi ved løse to ligninger efter $C_1$ og $C_2$.
%% Cell type:code id: tags:
``` python
x_fuld = C1 * E ** (2*t) + C2 * E ** t - 3/20 * sin(2*t) - 1/20 * cos(2*t)
solve([Eq(x_fuld.subs(t,0),0),Eq(diff(x_fuld,t).subs(t,0),1)])
```
%% Cell type:code id: tags:
``` python
plot(x_fuld.subs([(C1,1.25),(C2,-1.2)]),xlim = (-1,1), ylim = (-1,2))
```
%% Cell type:markdown id: tags:
## 2. Med dsolve
Nu finder vi samme løsning med $\text{dsolve}$
%% Cell type:code id: tags:
``` python
x = Function("x")
res = dsolve(Eq(diff(x(t),t,t) - 3*diff(x(t),t) + 2 * x(t),cos(2*t)),ics = {x(0) : 0, diff(x(t),t).subs(t,0) : 1})
res
```
%% Cell type:code id: tags:
``` python
plot(res.rhs,xlim = (-1,1),ylim = (-1,2))
```
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
%% Cell type:markdown id:a16f919a tags:
# Prikprodukt og krydsprodukt
Demo af Christian Mikkelstrup og Hans Henrik Hermansen
%% Cell type:code id:02fbdd15 tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id:3c6bf518 tags:
Prikprodukter og krydsprodukter bruges meget, især når forårssemesteret nåes. Heldigvis er det gode funktioner for begge dele indbygget i Sympy.
%% Cell type:markdown id:57df8306 tags:
## Dot og Cross
%% Cell type:markdown id:21f7179b tags:
I planen er giver 2 vektorer, find **prikproduktet** af dem:
%% Cell type:code id:23b1fa09 tags:
``` python
u = Matrix([-7,9])
v = Matrix([13,8])
u.dot(v)
```
%% Cell type:markdown id:908f8666 tags:
Eller generelt:
%% Cell type:code id:af28b92c tags:
``` python
a,b,c,d = symbols(['a', 'b', 'c', 'd'], real = True)
u = Matrix([a,b])
v = Matrix([c,d])
u.dot(v)
```
%% Cell type:markdown id:7468adb6 tags:
Med vektorer i rummet, kan vi også finde deres prikprodukt på samme måde
%% Cell type:code id:faa6eb05 tags:
``` python
u = Matrix([-7,9,2])
v = Matrix([13,8,-1])
u.dot(v)
```
%% Cell type:markdown id:5d99a5f5 tags:
Derudover kan vi også finde deres **krydsprodukt** ved
%% Cell type:code id:e28a9160 tags:
``` python
j = u.cross(v)
j
```
%% Cell type:markdown id:301e0d5a tags:
*Bemærk*: Prikproduktet er et **tal** (en skalar). Krydsproduktet giver en ny **vektor**, der er vinkelret på de to givne vektorer! Dette eftertjekkes ved
%% Cell type:code id:62b183b7 tags:
``` python
j.dot(u),j.dot(v)
```
%% Cell type:markdown id:33f98d7c tags:
Da prikproduktet giver 0, ses det at den nye vektor er vinkelret på begge vektorer i krydsproduktet. Som forventet.
%% Cell type:markdown id:1fd990e7 tags:
## Længder og vinkler
I dette eksempel udregner vi længden af og vinkler mellem vektorer ved brug af prikproduktet. Længden udregnes ved
\begin{equation}
|u|=\sqrt{u\cdot u},
\end{equation}
Bemærk at længden/normen $|u|$ oftest skrives $||u||_2$ i ingeniørvidenskaben. Længden af vektorene kan udregnes ved én af følgende 3 metoder:
%% Cell type:code id:11f9f90b tags:
``` python
u_1 = Matrix([1,2])
u_2 = Matrix([3,4])
```
%% Cell type:code id:50b5f60c tags:
``` python
# Gange sammen elementvis, sum sammen bagefter, kvadratroden til sidst
sqrt(sum(u_1.multiply_elementwise(u_1)))
```
%% Cell type:code id:53713bc1 tags:
``` python
# Indbyggede funktion for prikproduktet (dot product), og så kvadratroden
sqrt(u_2.dot(u_2))
```
%% Cell type:code id:eeecf6a4 tags:
``` python
# Indbyggede funktion, så vi kan bede om normen
# (2-normen er standard for vektorer, Frobenius normen for matricer)
u_2.norm() # er det samme som u_2.norm(2)
```
%% Cell type:markdown id:9158d5e7 tags:
Vinklerne mellem vektorene, $v$, udregnes udelukkende ved brug af ligheden fra *sætning 10.53*,
\begin{equation}
\cos(v)=\frac{u_1 \cdot u_2}{|u_1|\,|u_2|},
\end{equation}
hvor $u_1$ og $u_2$ er egentlige vektorer. Vi udregner vinkel mellem $u_1$ og $u_2$:
%% Cell type:code id:0bf3e613 tags:
``` python
#
cos_to_v = (u_1.dot(u_2))/(u_1.norm()*u_2.norm())
# arccos(x) fåes i python ved acos(x)
v = acos(cos_to_v)
v
```
%% Cell type:markdown id:2fd74371 tags:
hvilket cirka har værdien (målt i radianer)
%% Cell type:code id:abd3c7b3 tags:
``` python
N(v,5)
```
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
%% Cell type:markdown id: tags:
# Første ordens linære differentialligninger
Demo af Christian Mikkelstrup og Hans Henrik Hermansen
%% Cell type:code id: tags:
``` python
from sympy import *
init_printing()
```
%% Cell type:markdown id: tags:
Vi er givet følgende 1. ordens linære inhomogene differentialligning
\begin{equation*}
\frac{\text{d}}{\text{d}t}x(t) + x(t) = e ^ t
\end{equation*}
Vi ønsker at gøre følgende
1. Bestemme den fuldstændige løsning
2. Bestemme en løsning, der opfylder $x(0)=2$
3. Undersøge hvad startbetingelsen har af inflydelse på løsningen.
%% Cell type:markdown id: tags:
## Metode 1: Simuleret håndregning
Her vil anvende SymPy til at bruge panser-formlen på vores differntialligning
%% Cell type:code id: tags:
``` python
t,C = symbols('t C')
p = 1
P = integrate(p,t)
q = exp(t)
res = exp(-P) * (integrate(exp(P) * q,t) + C)
res
```
%% Cell type:markdown id: tags:
Vi har nu fundet den fuldstændige løsning. Nu skal vi finde en løsning, hvor $x(0)=2$. Da vi har fundet $x(t)$ kan vi indsætte 0 på $t$'s plads og løse for $C_1$
%% Cell type:code id: tags:
``` python
solve(Eq(res.subs(t,0), 2))
```
%% Cell type:markdown id: tags:
Altså har vi:
$x(t) = \frac{3}{2 e ^ t} + \frac{e ^ t}{2}$
Da dette er den eneste løsning er vi i overensstemmelse med *eksistens- og entydighedssætningen*.
Vi plotter løsningen:
%% Cell type:code id: tags:
``` python
res.subs(C,3/2).expand()
```
%% Cell type:code id: tags:
``` python
plot(res.subs(C,3/2),xlim=(-1,2),ylim=(0,5))
```
%% Cell type:markdown id: tags:
Vi kontrollerer at grafen faktisk går igennem punktet (0,2).
%% Cell type:markdown id: tags:
## Metode 2: Funktionen dsolve
Den næste metode gør brug af SymPy's indbyggede funktion $\text{dsolve}$, der løser differentialligninger.
Lige som vi definerer vores variable via funktionen $\text{symbols}$, når vi skal løse f.eks. ligninger,
så skal også definere funktionsvariabler, som repræsenterer ukendte funktioner, når vi løser differentialligninger.
Dette gøres via funktionen $\text{Function()}$,
%% Cell type:code id: tags:
``` python
x = Function('x')
dsolve(Eq(diff(x(t),t)+x(t),exp(t)))
```
%% Cell type:markdown id: tags:
Nu har vi den fuldstændige løsning. Vi kan også bruge SymPy til at finde en specifik løsning.
I det tilfælde skal blot bruge funktionsparameteren $\text{ics()}$, som er "initial conditions".
%% Cell type:code id: tags:
``` python
res = dsolve(Eq(diff(x(t),t) + x(t),exp(t)),ics = {x(0) : 2})
res
```
%% Cell type:markdown id: tags:
Her kan vi se, at vi får den samme løsning. Mon ikke også den ser ens ud, når vi plotter den.
%% Cell type:code id: tags:
``` python
plot(res.rhs,xlim=(-1,2),ylim=(0,5))
```
%% Cell type:markdown id: tags:
Som forventet!
Vi kan også bruge Pythons loops til at vise mange forskellige startbetingelser:
%% Cell type:code id: tags:
``` python
p1 = plot(show=False)
for x0 in range(11):
# Ny ligning med nye initial conditions
neweq = dsolve(Eq(diff(x(t),t) + x(t),exp(t)),ics = {x(0) : x0})
# Nyt plot, hvor vi plotter for t mellem -1 og 2
newplot = plot(neweq.rhs, (t, -1, 2), show=False)
# tilføj til vores plot
p1.extend(newplot)
p1.xlim = (-1,2)
p1.show()
```
%% Cell type:code id: tags:
``` python
```
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment