5.13. SymPy og Differentialligninger#

I denne Notebook vil vi kigge på brugen af SymPy til at løse differentialligninger. Formålet er for det første at introducere grundlæggende terminologi som både kan være relevant i LinALys og det følgende kursus MatF, hvor differentialligninger spiller en stor rolle. Det andet formål er at indføre nogle simple teknikker til at tjekke resultater opnået ved regning i hånden.

Til at starte med kører vi vores standard start-blok:

import sympy as sp                         # Importer sympy
from sympy.abc import x, y, a              # Vi definerer vores symbolske variabel og konstanter

5.13.1. Opskrivning og løsning af en simpel differentialligning#

Hvis vi er givet en differentialligning kan vi løse den næsten som en almindelig ligning i SymPy. Det er dog vigtigt at fortælle SymPy, at vores symboler nu repræsenterer en funktion og ikke en variabel. Dette gøres med sp.Function(funktionsnavn), hvor vi som navn oftest blot anvender betegnelserne “f”, “g” eller lignende. Vi fortæller SymPy, at eksempelvis f er en funktion af x ved at skrive f(x), og kan nu differentiere den ved brug af f(x).diff(x), som vi har også har gjort tidligere. For at formulere selve differentialligningen vil vi bruge sp.Eq til at angive at to udtryk er lig hinanden:

# Vi definerer en funktion:
f = sp.Function("f")

# Vi kan nu lave en differentialligning
DiffLigning = sp.Eq(f(x).diff(x), - a * f(x)) 
DiffLigning
\[\displaystyle \frac{d}{d x} f{\left(x \right)} = - a f{\left(x \right)}\]

Vi benytter nu sp.dsolve(differentialligning, funktionsnavn) til at bestemme en løsning til differentialligningen. Bemærk at det andet argument “funktionsnavn” er dvendigt for at angive hvad SymPy skal løse for, som med notationen ovenfor er f(x).

# Benyt nu dsolve
sp.dsolve(DiffLigning, f(x))
\[\displaystyle f{\left(x \right)} = C_{1} e^{- a x}\]

Vi får her den generelle løsning, idet konstanten \(C_1\) er ubestemt (mens \(a\) optræder i den oprindelige differentialligning. For at bestemme en løsning, der opfylder en given startbetingelse, bruges keywordet ics (som er en forkortelse for intial conditions). Formatet er et såkaldt dictionary, hvilket vil sige at startbetingelsen \(f(x_0) = c\) formuleres som {f(x_0): c}. Den samme differentialligning som ovenfor løses sammen med startbetingelsen \(f(0) = 5\) således :

sp.dsolve(DiffLigning, f(x), ics = {f(0): 5})
\[\displaystyle f{\left(x \right)} = 5 e^{- a x}\]

5.13.2. Andenordens differentialligninger#

Vi kan også løse mere avancerede differentialligninger, f.eks. en andenordens differentialligning. Her løses følgende:

DiffLigning = sp.Eq(f(x).diff(x, 2) - a * f(x).diff(x), - a ** 2 * f(x))
DiffLigning
\[\displaystyle - a \frac{d}{d x} f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = - a^{2} f{\left(x \right)}\]
sp.dsolve(DiffLigning, f(x))
\[\displaystyle f{\left(x \right)} = C_{1} e^{\frac{a x \left(1 - \sqrt{3} i\right)}{2}} + C_{2} e^{\frac{a x \left(1 + \sqrt{3} i\right)}{2}}\]

Skal vi finde en partikulær løsning til en andenordens differentialligning, skal vi benytte to stykker information til at bestemme \(C_1\) og \(C_2\), f.eks. løsningens værdi i to punkter eller en kombination af et punkt og løsningkurvens hældning i punktet.

Skal man angive en startbetingelse som f.eks. \(f'(0) = 3\), da man skal differentiere i forhold til x og så indsætte \(x=0\). Dette gøres ved f(x).diff(x).subs(x, 0). Vi kan samle flere startbetingelser ved blot at seperere dem med et komma. Lad os løse overstående differential ligning med startbetingelserne: \(f'(0) = 3\) og \(f(0) = 10\).

sp.dsolve(DiffLigning, f(x), ics = {f(0): 10, f(x).diff(x).subs(x, 0): 3})
\[\displaystyle f{\left(x \right)} = \left(5 - \frac{5 \sqrt{3} i}{3} + \frac{\sqrt{3} i}{a}\right) e^{\frac{a x \left(1 - \sqrt{3} i\right)}{2}} + \left(5 + \frac{5 \sqrt{3} i}{3} - \frac{\sqrt{3} i}{a}\right) e^{\frac{a x \left(1 + \sqrt{3} i\right)}{2}}\]

5.13.3. Check solution#

Vi vil her introducere hvordan man kan tjekke, om vi har fundet en korrekt løsning til en differentialligning. Her benytter vi sp.checkodesol (som er en sammentrækning af check ordinary differential equation solution). Vi skal angive en differentialligning (f.eks. formuleret ved hjælp af sp.Eq ligesom ovenfor), og den løsning, som vi enten har fundet ved hændregning eller eventuelt gættet. Funktionen returnerer True hvis løsningen er rigtig, eller False sammen med de led, som var tilbage ved indsættelse af løsningen i differentialligningen. Lad os tjekke, om vi kan gætte en løsning til den differentialligning, der beskriver den harmoniske oscillator.

HamOsc = sp.Eq(f(x).diff(x, 2), -a * f(x))
HamOsc
\[\displaystyle \frac{d^{2}}{d x^{2}} f{\left(x \right)} = - a f{\left(x \right)}\]

Vi gætter nu på at løsningen kunne være en cosinus-funktion:

sp.checkodesol(HamOsc, sp.cos(x))

(False, (1 - a)*cos(x))

Vi ser at der i hvert tilfælde er et problem med konstanten \(a\), hvilket motiverer os til at opdatere vores gæt:

sp.checkodesol(HamOsc, sp.cos(sp.sqrt(a) * x))

(True, 0)

Selv hvis man ikke kan gætte løsninger, kan metoden være en god hjælp til den skriftlige eksamen, hvor Python/SymPy er et tilladt hjælpemiddel. Dette gælder især de opgaver, hvor man skal finde en løsning og regne videre med resultatet i næste delspørgsmål. Her kan det være en god investering at checke løsningen inden man går videre.