Ekstra opgaver

Den simpleste måde at løse differentialligninger på blev introduceret af den schweiziske matematiker Leonhard Euler (1707-1783) og kaldes derfor Euler-metoden. Se på differentialligningen:

Kendes en værdi af den afhængige variabel, $y(t_0)=y_0$, til en bestemt værdi af den uafhængige variabel, "tiden" $t_0$, kan senere værdier bestemmes med Euler-metoden, der iterativt frembringer løsningen ved at tage små tidsskridt fra det $n$'te til det $(n+1)$'te punkt på en diskret repræsentation af kurven:

I Euler-metoden beregnes den næste funktionsværdi alene på basis af den aktuelle differentialkvotient. Dvs. tid og funktionsværdi fremskrives til første orden i en skridtlængde, $\Delta t$. Denne metode er eksakt hvis $f$ er uafhængig af $t$ og $y$, men hvis $f$ varierer, akkumulerer Euler-metoden hurtigt fejl. Disse kan reduceres ved anvendelse af højereordensdifferentialigningsløsere.

En måde at forbedre Euler-metoden på er vist i figuren ovenfor. Her beregnes først én hældning, $s_1=dy/dt$, for det aktuelle punkt, dernæst en ny hældning, $s_2$, "et halvt skridt fremme" for så endeligt at tage hele skridtet baseret på denne, nye hældning:

En sådan fremgangsmåde kaldes en Runge-Kutta-metode efter dens opfindere, C. Runge og M.W. Kutta. Metoden er af 2. orden da den akkumulerede fejl vokser proportionalt med $\Delta t^2$. Betragt differentialligningen,

for $0\le t\le 5\pi$. a) Løs først differentialligningen efter Euler-metoden:

implementeret som en for-løkke, noget i stil med dette ufuldstændige program, hvor t og y beskriver $(t_0,y_0)$ til at starte med og så løbende opdateres:

ts = [t]
ys = [y]

for _ in np.arange(t,tmax,delta_t):
   y += ...
   t += ...
   ts.append(t)
   ys.append(y)

Plot resultatet for 3 forskellige værdier af $\Delta t$, fx $0.01$, $0.05$ og $0.25$.

b) Gentag spgm. a) idet du løser differentialligningen efter en 2. ordens Runge-Kutta metode:

c) Løs nu differentialligningen under anvendelse af solve_ivp, hvor parameteren, der tilsvarer $\Delta t$ hedder max_step.

d) Skriv til sidst din løsning til spgm. b) som en funktion, my_ivp_solver, der kaldes med tre positionelle argumenter og et keyword-argument:

1. en funktionshenvisning til en funktion, $f(t,y)$, der returnerer $\displaystyle\frac{dy}{dt}$
2. en liste, $[t_{start}, t_{slut}]$, der angiver start og slutværdier for tiden, $t$
3. en liste med ét tal, $y_0$, der er startværdien for $y$, dvs. $y(t_{start})$.
4. et keyword-argument, max_step, der angiver værdien af $\Delta t$.

Lad endvidere din funktion aflevere de to lister med tider og funktionsværdier på samme måde som solve_ivp, hvilket kan gøres ved at sætte dem på som attributter til en ellers "tom" funktion:

    ...
        f = lambda x: None
        f.t = np.array(ts)
        f.y = np.array([ys])
    return f

(Dette er lidt et hack, hvor det udnyttes at man kan sætte attributter på funktioner. f kunne have været et andet python objekt end lige en funktione, men fx ikke en integer eller en liste.) Hvis du skriver koden rigtigt skal du kunne udskifte solve_ivp med my_ivp_solver i din kode til spørgsmål c og få samme plot som i delspørgsmål b.

Et ustrækkeligt tov med masse per længdeenhed, $\lambda$, er ophængt lige højt i begge ender. Under indvirkning af tyngdekraften vil det hænge som på figuren og følge det der betegnes kædelinien (Eng: catenary):

Kraften i tovet, snorkraften, virker altid langs tovets retning, hvilket betyder at snorkraften midt i tovet (dvs. for $(x,y)=(0,0)$, jf. figuren) peger i $\pm x$-retningen. Størrelsen af snorkraften betegnes $S_0$. Går man en længde $s$ ud ad tovet til et punkt $(x_s,y_s)$ vil al vægten af tovet fra midten til $s$ mærke en tyngdekraft, $s\lambda g$ i $-y$-retningen, hvor $g$ er tyngdeaccelerationen. Kraftligevægt for dette stykke af tovet betyder derfor at snorkraften i $s$ må være $S_0$ i $+x$-retningen og $s\lambda g$ i $+y$-retningen. Da snorkraften kun kan pege i tovets retning vides nu at tovets hældning som funktion af længden af tov, $s$, fra tovets midte opfylder:

hvor $a$ er introduceret for at forkorte dette og følgende udtryk. Vi har tidligere i kurset arbejdet med kurveintegraler og længden af en kurve:

Benyttes nu differentialregningens hovedsætning fås følgende udtryk for $\displaystyle\frac{ds}{dx}$:

Indsættes $\displaystyle\frac{dy}{dx}=\frac{s}{a}$, heri og uddrages det reciprokke står man med:

Med kæderegelen, $\displaystyle\frac{dy}{ds}=\frac{dy}{dx}\frac{dx}{ds}$, kan man nemt udlede et tilsvarende udtryk for $\displaystyle\frac{dy}{ds}$:

Vi har altså $x$- og $y$-koordinaterne for kædelinien beskrevet ved hver deres 1. ordensdifferentialligning, hvor den uafhængige variable er den samme, nemlig længden af kædelinien, $s$. Vi har dermed to muligheder for at finde og plotte kædelinien:

1. Vi kan løse de to 1. ordensdifferentialligninger hver for sig med to kald af solve_ivp og få løsningerne for hhv $x(s)$ og $y(s)$ for de samme værdier af $s$.
2. Vi kan løse de to 1. ordensdifferentialligninger samtidigt som om de var koblede differentialligninger (selvom de ikke er det) med ét kald af solve_ivp. Hermed bliver det automatisk opfyldt at $x(s)$ og $y(s)$ fås for de samme værdier af $s$.

Dette skal nu prøves i praksis, idet $a=2$ benyttes og en passende startbetingelse aflæses fra figuren. a) Benyt den første metode og tjek at $s$-værdierne er de samme. Benyt evt t_eval-argumentet i kaldet af solve_ivp for selv at styre hvilke $s$-værdier, der benyttes. Plot kædelinien. b) Benyt den anden metode og plot oven i svaret fra spørgsmål a). c) Plot den analytiske løsning, $\displaystyle y=a\cosh(\frac{x}{a})-a$ for $a=2$, og konstatér at den også ligger oven i de numeriske løsninger.