Numerisk løsning af differentialligninger

I størstedelen af de fysiske problemer man løser med simulering, indgår numerisk integration. Det er naturligt da næsten alle fysiske modeller er formuleret som differentialligninger. Den simpleste form for numerisk integration kaldes "Euler metoden", og er nem at forstå, men også meget upræcis, da den kun udnytter ændringen i begyndelsen af et skridt. Til seriøse beregninger er metoden uegnet, men den er meget illustrativ da hvert skridt er ekstrapoleret lineært fra det forrige.
En noget bedre teknik, som ofte er brugt er Runge-Kutta metoden. Runge-Kutta (RK) tager som Euler metoden et skridt af en numerisk størrelse fra den forrige værdi. Men i stedet for kun at udnytte ændringen i begyndelsen af skridtet, estimerer den ændringen ud fra et eller flere "prøve punkter" i skridt-intervallet, som vist her:

Euler - Næste trin afhænger kun af forrige.

Runge-Kutta - Næste trin afhænger også af værdien mellem den forrige og næste værdi.

Euler metoden

Som en del af de eksperimentielle øvelser i Mek1 og Mek2 har vi simuleret forskellige systemer ved at løse Newtons 2. lov numerisk:

\begin{equation} \label{eq:n2} \vec{F} = m \frac{d^2 \vec{r}}{d t^2} \end{equation}

Teknikken vi har benyttet hedder `Eulers metode' og bygger på at man for et lille interval $h$ approksimerer ændringen i en variabel \(y_{n+1}\) på baggrund af den tidligere værdi \(y_n\) og funktionsudtrykket for den resterende del af ligningen \(f(t, y)\).

\begin{align} \label{eq:euler} y_{n+1} = y_n + h \ f(t_n, y_n) \end{align}

Boksen herunder giver et eksempel på hvordan man løser en differentialligning med Eulers metode. Eulers metode virker kun for 1-ordens differentialligninger. Ligning~\ref{eq:n2} er en 2-ordens differentialligning, som kan omskrives til et sæt af 1-ordensligninger ved substitution,

\begin{align} m \frac{d v}{d t} &= F \\ \frac{d r}{d t} &= v \end{align}

Som vi løser for hhv. \(d v\) og \(d r\),

\begin{align} d v &= F/m \ dt\\ d r &= v \ dt \end{align}

Til sidst erstatter vi \(d t\) med en endelig skridt størrelse \(h\) samt adderer ændringen til det tidligere tidsskridt og får de to ligninger som beskriver udviklingen i hastighed \(v\) og position \(r\) for et objekt med massen \(m\) udsat for kraften \(F\),

\begin{align} \label{eq:n2_euler} r_{n+1} = r_n + h \ v_n \\ v_{n+1} = v_n + h \ F/m. \end{align}

Bemærk at \(r_{n+1}\) opdateres af \(v_n\) og ikke \(v_{n+1}\), det kan være fristende at bytte om på rækkefølgen af de to ligninger i sin kode, men så er hastighedsestimatet altså \(h\) sekunder foran og positionen og hastigheden vil blive forskudt en smule.

Eksempel på Eulers metode

Løs lig.~\ref{eq:euler_ex} med Eulers metode. \begin{equation} \label{eq:euler_ex} \frac{dy}{dt} = y \quad y(0) = 1 \end{equation} Først omskrives udtrykket således at vi isolerer \(d y\) som er ændringen vi er interesseret i, \begin{equation} d y = y \ d t \end{equation} For at vise at vi nu arbejder med endelige tal erstatter vi \(d t\) med \(h\). I forhold til lig.~\ref{eq:euler} er \(f(t_n, y_n) = y_n\), her er ingen tidsafhængighed, men det næste skridt \(y_{n+1}\) afhænger af det forrige, \begin{equation} y_{n+1} = y_n + y_n \ h \quad y_0 = 1. \end{equation} Begyndelsesbetingelsen \(y(0)=1\) betyder at vi starter med at sætte vores variabel \(y=1\), hvis vi sætter intervalstørrelsen \(h=0.1\) kan vi beregne de første \(10\) værdier: \begin{align} y_0 &= 1 \notag \\ y_1 &= y_0 + y_0 \ h = 1 + 1 \times 0.1 = 1.1 \notag \\ y_2 &= y_1 + y_1 \ h = 1.1 + 1.1 \times 0.1 = 1.21\notag \\ y_3 &= y_2 + y_2 \ h = 1.21 + 1.21 \times 0.1 = 1.331\notag \\ &\vdots \notag \\ y_{10} &= y_9 + y_9 \ h = 2.36 + 2.36 \times 0.1 = 2.594 \end{align} I MATLAB vil ovenstående se således ud,
	y(1) = 1; % begyndelesbetingelse (liste)
	h = 0.1; % interval
	for i=1:10 % antal skridt i alt
	    y(i+1) = y(i) + h * y(i); % beregn næste skridt
	end
Efter løkken er gennemløbet vil y være en liste med \(11\) værdier, som hver svarer til \(x(n)=0...h \times n\) Formuleret i MATLAB, med en \(x\) værdi i stedet for \(n\),
y(1) = 1; % begyndelesbetingelse (liste)
x(1) = 0; % begyndelsespunkt
x_end = 1; % endepunkt
h = 0.1; % interval størrelse
n = 1; % tæller
while x(end) < x_end
    x(n+1) = x(n) + h; % opdater position på x-aksen
    y(n+1) = y(n) + h * y(n); % beregn næste skridt
    n = n + 1; % opdater tæller
end

plot(x, y,'r') % plot det numeriske resultat
hold on
% Analytisk løsning med dsolve()
y_a = dsolve('Dy=y','y(0)=1','x'); % Løsning: @(x)exp(x)
ezplot(y_a,[x(1), x_end]) % plot det analytiske resultat

Figuerne herunder viser resultatet af koden i listing~\ref{lst:eulerfull} for to forskellige \(h\)-værdier: \(h=0.2, 0.01\). Det er tydeligt at jo mindre skridtstørrelser man vælger jo bedre en approksimation får man af den analytiske løsning.

Opgaver i Euler metoden

  1. Udfør programmet i listing~\ref{lst:eulerfull} men for en anden ligning, f.eks. \(dy/dx=y^2, y(0)=-5\).

Runge-Kutta metoden

Euler metodens store fordel er at den er gennemskuelig. Hver iteration svarer til et tidsskridt, det svarer til at man spoler Newtons kalkulus `baglæns' og i stedet for at tage uendeligt små skridt, så gør man dem endeligt store. Desværre er metodens præcision meget begrænset da antallet af beregninger hurtigt bliver uoverkommeligt selv med relativt store skridtstørrelser. Gennem tiderne er der udviklet andre teknikker som på forskellig vis forbedrer præcisionen uden at det kræver væsentligt flere beregninger. En af disse metoder kaldes Runge-Kutta. Ideen bag denne teknik er at man i stedet for at antage at en funktion ikke ændrer sig næveværdigt mellem \(y(h)\) og \(y(2 \times h)\) laver \(4\) estimater i området og benytter et vægtet gennemsnit af dem som det endelige bud på hvad \(y(2 \times h)\) er.

\begin{align} \label{eq:rk4} k_1 &= h f(t_n, y_n) \quad \text{`Eulers metode'}\\ k_2 &= h f(t_n + \frac{1}{2} h, y_n + \frac{1}{2} k_1)\\ k_3 &= h f(t_n + \frac{1}{2} h, y_n + \frac{1}{2} k_2)\\ k_4 &= h f(t_n + h, y_n + k_3)\\ \\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4) \quad \text{`Runge-Kutta $\mathcal{O}(4)$'}\\ t_{n+1} &= t_n + h \end{align}

Eksempel på Runge-Kutta metoden

Løs lig.~\ref{eq:rk_ex} med Runge-Kutta metoden til 4. orden.

\begin{equation} \label{eq:rk_ex} \frac{dy}{dt} = y \quad y(0) = 1 \end{equation}

Først omskrives udtrykket således at vi isolerer \(dy\) som er ændringen vi er interesseret i,

\begin{equation} d y = y \ d t \end{equation}

For at vise at vi nu arbejder med endelige tal erstatter vi \(d t\) med \(h\). I forhold til lig.~\ref{eq:rk_ex} er \(f(t_n, y_n) = y_n\), her er ingen tidsafhængighed, men det næste skridt \(y_{n+1}\) afhænger af det forrige,

\begin{align} k_1 &= h y_n\\ k_2 &= h (y_n + \frac{1}{2} k_1)\\ k_3 &= h (y_n + \frac{1}{2} k_2)\\ k_4 &= h (y_n + k_3)\\ \\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4) \quad y_0 = 1 \\ \end{align}

Begyndelsesbetingelsen \(y(0)=1\) betyder at vi starter med at sætte vores variabel \(y=1\), hvis vi sætter intervalstørrelsen \(h=0.1\) kan vi beregne den første værdi:

\begin{align} k_1 &= h y_n = 0.1 \times 1 = 0.1\\ k_2 &= h (y_n + \frac{1}{2} k_1) = 0.1 \times (1 + \frac{1}{2} \ 0.1) = 0.105\\ k_3 &= h (y_n + \frac{1}{2} k_2) = 0.1 \times (1 + \frac{1}{2} \ 0.105) = 0.10525\\ k_4 &= h (y_n + k_3) = 0.1 \times (1 + 0.10525) = 0.110525\\ \\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4) = 1.1051708 \end{align}

I MATLAB vil ovenstående se således ud,

y_rk(1) = 1; % runge kutta begyndelse
h = 0.5; % interval størrelse
for n=1:10
    k1 = h * y_rk(n);
    k2 = h * (y_rk(n) + 1/2*k1);
    k3 = h * (y_rk(n) + 1/2*k2);
    k4 = h * (y_rk(n) + k3);
    y_rk(n+1) = y_rk(n) + 1/6 * (k1 + 2*k2 + 2*k3 + k4); % RK4
end

I eksemplet herunder sammenligner vi Runge-Kutta, Euler metoden og den analytiske løsning,

y(1) = 1; % euler begyndelse
y_rk(1) = 1; % runge kutta begyndelse
x(1) = 0; % begyndelsespunkt
x_end = 4; % endepunkt
h = 0.5; % interval størrelse
n = 1; % tæller
while x(end) < x_end
    x(n+1) = x(n) + h; % opdater position på x-aksen

    k1 = h * y_rk(n);
    k2 = h * (y_rk(n) + 1/2*k1);
    k3 = h * (y_rk(n) + 1/2*k2);
    k4 = h * (y_rk(n) + k3);

    y(n+1) = y(n) + k1; % Euler
    y_rk(n+1) = y_rk(n) + 1/6 * (k1 + 2*k2 + 2*k3 + k4); % RK4

    n = n + 1; % opdater tæller
end
h0 = ezplot(dsolve('Dy=y','y(0)=1','x'),[x(1), x_end]) % plot det analytiske resultat
hold on
h1 = plot(x, y,'r-o') % plot det numeriske resultat
h2 =plot(x, y_rk,'k-o') % plot det numeriske resultat
legend([h1 h2 h0],['Euler metoden h=' num2str(h)],['Runge-Kutta h=' num2str(h)],'Analytisk f(x)=e^{x}')

Som det fremgår af plottet er RK4 metoden langt mere præcis i forhold til Euler metoden ved samme skridt-størrelse (\(h=0.5\)) for den valgte funktion. Til simuleringer som skal køre længe og derfor tage mange integrationstrin har det stor betydning at estimatet er præcist, da fejlen forværres kommulativt.

Opgaver i Runge-Kutta integration

  1. Implementer RK4 metoden som vist ovenfor. Plot forskellige \(h\) værdier.
  2. Lav et plot der viser den relative fejl mellem de tre løsninger \(fejl = \frac{y_{ana} - y_{rk}}{y_{ana}}\) som funktion af x.
  3. På fejl-plottet kan man se hvordan afvigelsen fra den analytiske værdi vokser med en hvis størrelse per punkt, giv en kvalitativ forklaring på hvor fejl opstår i de numeriske metoder.
  4. Hvad skal skridtstørrelsen være for Euler metoden før resultatet svarer til RK4 med \(h=0.5\)?
  5. Implementer Newtons bevægelsesligninger med RK4 metoden (lig.~\ref{eq:n2_euler}), med følgende begyndelsesbetingelser: \(F=1/r\) [N], \(m=10\) [kg], \(r(0)=0\) [m], \(v(0)=10\) [m/s]. Hvor befinder legemet sig efter 10 sek? (Hint: Som med Eulers metode skal der være to sæt af \(k_{1,2,3,4}\) et for hver differentialligning \(v(t)\) og \(r(t)\). Husk også rækkefølgen på opdateringen af \(v_{n+1}\) og \(r_{n+1}\)).