\documentclass[a4paper]{scrartcl}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{ngerman}
\usepackage[dvips]{graphicx}
\usepackage{ulem}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{array}
\usepackage{booktabs}
\usepackage{listings}
\usepackage{color} % f{\"u}r den grauen Hintergrund


%\usepackage{endnotes}

\setlength{\parindent}{0 em}

\definecolor{Gray}{gray}{0.9}

\lstset{basicstyle=\ttfamily\scriptsize,
backgroundcolor=\color{Gray}, numbers=left, numberstyle=\tiny,
stepnumber=1, numbersep=5pt}
%\usepackage{endnotes}

\setlength{\parindent}{0 em}


\DeclareGraphicsExtensions{.eps}

%\DeclareGraphicsRule{.gif}{eps}{}{`convert #1 'eps:-'}

\begin{document}
\begin{center}
\textbf{\huge{Runde 3, Beispiel 19}}

\Large{LVA 118.181, {\"U}bungsrunde 3, 03.11.}

\large{Markus Nemetz, \emph{markus.nemetz@tuwien.ac.at}, TU Wien,
01.11.2006}
\end{center}

\section{Angabe}
Man vergleiche das Euler-Verfahren, das Verfahren von Heun und das
Verfahren von Runge- Kutta f{\"u}r das AWP $y = 1 + \frac{y}{x},
y(1)=1$ an der Stelle $x=1.6$, wobei $n = 3$ gew{\"a}hlt werden soll.
L{\"o}se das AWP auch exakt und ermittle jeweils den relativen Fehler
f{\"u}r die einzelnen Verfahren.

Anmerkung: Nach M{\"o}glichkeit programmiere man die einzelnen
Verfahren. Wenn 'mit der Hand' gerechnet werden muss, w{\"a}hle $n =
1$.

\section{Theoretische Grundlagen: Euler-Verfahren}
Das Euler-Verfahren ist ein numerisches Vefahren zur L{\"o}sung von
Anfangswertproblemen. Die Idee ist, denn Differentialquotienten
$y' = \frac{\partial y}{\partial x}$ durch den
Differenzenquotienten $\frac{\triangle y}{\triangle x} = \frac{y(x
+ h) - y(x)}{h}$ zu ersetzen ($h$ ist die Schrittweite):
\begin{gather*}
    f(x,y(x)) = y'(x) = \frac{\partial y}{\partial x} \approx \frac{\triangle y}{\triangle x} = \frac{y(x
+ h) - y(x)}{h}
\end{gather*}
Die Aufl{\"o}sung nach $y(x+h)$ ergibt
\begin{gather*}
    y(x + h) \approx y(x) + hf(x,y(x))
\end{gather*}
Wenn wir nun den bekannten Anfangswerk $y_0 = y(x_0)$, eine kleine
Schrittweite $h$ w{\"a}hlen und $x=x_0$ und $y=y_0$ setzen, erhalten
wir mit $y_1$ eine N{\"a}herung f{\"u}r den exakten Wert $y(x_0 +h)$:
\begin{gather*}
    y(x_1) = y(x_0 + h) \approx y_0 + fh(x_0,y_0) =: y_1
\end{gather*}
Daraus k{\"o}nnen wir $y_2$ errrechnen usw. worazs such folgende
allgemeine Formel ergibt:
\begin{gather*}
    y_{n+1} = y_n + hf(x_n,y_n)
\end{gather*}

\section{Theoretische Grundlagen: Heun-Verfahren}
Das Heun-Verfahren, benannt nach Karl Heun, ist ein einfaches
Verfahren zur numerischen L{\"o}sung von Anfangswertproblemen. Es ist
ein Einschrittverfahren und geh{\"o}rt zu der Klasse der
Runge-Kutta-Verfahren.

Im Gegensatz zum Explizites Euler-Verfahren erfolgt die N{\"a}herung
{\"u}ber ein Trapez und nicht {\"u}ber ein Rechteck.

Zur numerischen L{\"o}sung des Anfangswert-Problems:
\begin{gather*}
    \dot{x}=f(t,x), \quad \quad x(t_0)=x_0
\end{gather*}
f{\"u}r eine gew{\"o}hnliche Differentialgleichung mit dem Verfahren von
Heun w{\"a}hle man eine Diskretisierungs-Schrittweite $h>0$, betrachte
die diskreten Zeitpunkte
\begin{gather*}
    t_k=t_0+kh, \quad \quad k=1,2,\dots </math>
\end{gather*}
und berechne zun{\"a}chst analog zum expliziten Euler-Verfahren
\begin{gather*}
    x^{[P]}_{k+1}=x_k+hf(t_k,x_k) \quad,\quad k=1,2,\dots
\end{gather*}
und dann
\begin{gather*}
    x_{k+1}=x_k+\frac{1}{2}h(f(t_k,x_k)+f(t_{k+1},x^{[P]}_{k+1}))\quad,\quad k=1,2,\dots
\end{gather*}
Die $x_i$ sind die N{\"a}herungswerte der tats{\"a}chlichen
L{\"o}sungsfunktion $x(t)$ zu den Zeitpunkten $t_i$.

$h$ bezeichnet man als Schrittweite; Verkleinert man die
Schrittweite so wird der Verfahrensfehler kleiner (sprich die
$x_i$ liegen n{\"a}her am tats{\"a}chlichen Funktionswert $x(t_i)$). Der
globale Fehler des Verfahren von Heun geht mit $h^2$ gegen Null
man spricht auch von Konvergenzordnung 2.

\section{Theoretische Grundlagen: Runge-Kutta-Verfahren}
klassische Runge-Kutta-Verfahren  ist ein explizites 4-stufiges
Einschrittverfahren der numerischen Mathematik zur n{\"a}herungsweisen
L{\"o}sung von Anfangswertproblemen. Der Name stammt von den deutschen
Mathematikern Carl Runge und Martin Wilhelm Kutta-

Das klassische Runge-Kutta-Verfahren verwendet - wie die weitaus
meisten numerischen L{\"o}sungsverfahren f{\"u}r Differentialgleichungen -
den Ansatz, Ableitungen (Differentialquotient) durch (endliche)
Differenzenquotienten zu ersetzen. Die dabei bei nichtlinearen
Funktionen notwendigerweise auftretenden Fehler (es werden
s{\"a}mtliche h{\"o}heren Glieder der Taylor-Entwicklung vernachl{\"a}ssigt)
k{\"o}nnen durch geeignete Kombinationen verschiedener
Differenzquotienten teilweise kompensiert werden. Das
Runge-Kutta-Verfahren ist nun eine solche Kombination, die
Diskretisierungsfehler bis zur dritten Ableitung kompensiert.

Es sei
\begin{gather*}
    y'=f(x,y)
\end{gather*}

eine (gew{\"o}hnliche) Differentialgleichung 1.Ordnung mit den
Anfangsbedingungen
\begin{gather*}
    y(x_{0})=y_{0}
\end{gather*}
und sei weiter $h$ die gew{\"u}nschte Schrittweite.

Dann wird der (angen{\"a}herte) Wert von $y(x_0+h)=y(x_1)$ nach Runge
wie folgt errechnet:

        \begin{center}
            \includegraphics[scale=0.6]{runge-tab.eps}
        \end{center}


Hinweis: In der Literatur wird h{\"a}ufig an Stelle der $y'_X$ direkt
die Schrittdifferenz $k_X=h.y'_X$ angegeben. Das ist insbesondere
hinsichtlich der in den folgenden Abschnitten angegebenen
Anwendung auf Differentialgleichungen h{\"o}herer Ordnung unpraktisch
und wird darum (weil es die schematische Rechenarbeit unterbricht
und auch sonst keinen Gewinn bringt) hier nicht verwendet.
\begin{gather*}
\begin{matrix}
y'_0 &=& f(x_0,y_0)\\
y_A &=& y_0+\frac{h}{2} \cdot y'_0\\
y'_A &=& f(x_0+\frac{h}{2},y_A)\\
y_B &=& y_0+\frac{h}{2} \cdot y'_A\\
y'_B &=& f(x_0+\frac{h}{2},y_B)\\
y_C &=& y_0+h \cdot y'_B\\
y'_C &=& f(x_0+h,y_C)\\
y_1 &=& y_0+h \cdot \frac{y'_0+2(y'_A+y'_B)+y'_C}{6}\\
\end{matrix}
\end{gather*}
Mit den neuen Werten $x_1$ und $y_1$ kann dann der n{\"a}chste
Rungeschritt durchgef{\"u}hrt werden. Die erzielte Genauigkeit liegt -
bei gen{\"u}gend glattem $f(x,y)$ - in der Gr{\"o}{\ss}enordnung von $h^4$.



\section{L{\"o}sung des Beispiels}
\begin{gather*}
    y = 1 + \frac{y}{x} \qquad \Rightarrow \qquad y - \frac{y}{x}
    = 1
\end{gather*}
Eine M{\"o}glichkeit zur L{\"o}sung besteht in der \textbf{Verwendung
einer Formel}: Lineare inhomogene Differentialgleichungen 1.
Ordnung in der Form $y' + p(x)y = r$ ($r$ ist St{\"o}rfunktion,
singul{\"a}r oder nur von $x$ abh{\"a}ngig) k{\"o}nnen mit folgender Formel
aufgel{\"o}st werden:
\begin{gather*}
    h = \int p(x) \, dx \\ y(x) = e^{-h} (\int e^hr \, dx + c)
\end{gather*}
(Quelle: Kreyszig, Advanced Engineering Mathematics, 9.Aufl.,
S.26ff.)
\begin{gather*}
    r=1, \qquad p = -\frac{1}{x}, \qquad h  \int p \, \,\, dx = -\ln
    x\\
    y(x) = \underbrace{e^{\ln x}}_{x} \cdot (\int \underbrace{e^{-\ln x}}_{\frac{1}{x}} \cdot 1 + c) \\
    y(x) = x \cdot (\int \frac{1}{x} \cdot 1 + c) \\
    y(x) = x \cdot (\ln x + c) \\
    y(x) = x\cdot \ln x + c
\end{gather*}
Spezielle L{\"o}sung f{\"u}r $y(1)=1$:
\begin{gather*}
    y(x) = x\cdot \ln x + x
\end{gather*}

Ausprogrammiert in MATLAB:

\begin{lstlisting}[caption={Definition der Funktion}]
%  Funktion phi=f(t,y) ist die rechte Seite von DGL y'(t)=f(t,y).
%  Laut der MATLAB-Syntax muss dann die Datei den Name "f.m" haben.

  function phi=f(t,y)

    phi=1 + t./y;
\end{lstlisting}

\begin{lstlisting}[caption={numdgl.m}]
%  Numerische Loesung der DGL y'(t)=f(t,y) mit dem
%    *  expliziten Euler-Verfahren
%    *  verbesserten Euler-Verfahren (Euler-Heun)
%    *  Runge-Kutta-Verfahren 4. Ordnung.
%
%  Die Funktion f(t,y) ist in der datei "f.m" angegeben.

   clear all;

   h=0.1;    %  Schrittweite
   N=10;    %  Anzahl der Schritte
   t=1.6;        %  Anfangszeit
   T=[t:h:t+h*N];    %  Array fuer Graphen
   y0=1;    %  Anfangswert
   YEX=t.*log(t) +t;    %  Exakte Loesung


   yee=y0;    %  yee ist die numerische Loesung mit dem expliziten Euler-Verfahren
   YEE(1)=y0;    %  YEE ist das Array fuer Graphen

   yeh=y0;     %  yeh ist die numerische Loesung mit dem Euler-Heun-Verfahren
   YEH(1)=y0;    %  YEH ist das Array fuer Graphen

   yrk=y0;     %  yrk ist die numerische Loesung mit dem Runge-Kutta-Verfahren
   YRK(1)=y0;    %  YRK ist das Array fuer Graphen


   for n=1:N


%  Euler-Explizit

     yee=yee+h*f(t,yee);    %  Ein Schritt des expliziten Euler-Verfahrens
     YEE(n+1)=yee;

%  Euler-Heun

     k1=f(t,yeh);
     k2=f(t+h/2, yeh+h/2*k1);
     yeh=yeh+h*k2;        %  Ein Schritt des Euler-Heun-Verfahrens

     YEH(n+1)=yeh;


%  Runge-Kutta

     k1=f(t,yrk);
     k2=f(t+h/2, yrk+h/2*k1);
     k3=f(t+h/2, yrk+h/2*k2);
     k4=f(t+h, yrk+h*k3);
     yrk=yrk+h/6*( k1+2*(k2+k3)+k4 );        %  Ein Schritt des Runge-Kutta-Verfahrens
     YRK(n+1)=yrk;



%  Naechster Zeitschritt

     t=t+h;

   end


%  Graphen

   plot(T, YEX, 'k');
   hold on
   plot(T,YEE,'r');
   plot(T,YEH,'g');
   plot(T,YRK,'b');
   hold off

   title('Numerische Loesung von DGL');
   xlabel('x');
   ylabel('Loesung y');
   legend('Exakte Loesung','Euler-Explizit','Euler-Heun','Runge-Kutta');
\end{lstlisting}

\end{document}
