\documentclass[10pt,a4paper]{article}

%% Polonizacja
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}
\usepackage{polski}
\usepackage{indentfirst}

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts}
\usepackage{latexsym}

%% listing (na razie nieużywane)
%\usepackage{listings}
%\lstloadlanguages{C}

%% Inne
\usepackage{comment}

%% Hiperłącza (powinien byc ostatnim ładowanym pakietem)
\usepackage[dvips]{hyperref}



%% -------------------------------------------------------------------------
%% Twierdzenia
\newtheorem{thm}{Twierdzenie}[section]
\newtheorem{df}{Definicja}[section]

\DeclareMathOperator{\sign}{sign}

\numberwithin{equation}{section}

\newcommand*{\n}[1][n]{$#1$\nobreakdash-\hspace{0pt}}
\newcommand*{\abs}[1]{\ensuremath{\lvert #1 \rvert}}
\newcommand*{\norm}[1]{\ensuremath{\lVert #1 \rVert}}
\newcommand*{\ABS}[1]{\ensuremath{\left\lvert #1 \right\rvert}}
\newcommand*{\NORM}[1]{\ensuremath{\left\lVert #1 \right\rVert}}
\newcommand*{\vect}[1]{\ensuremath{\mathbf{#1}}}
\newcommand*{\matr}[1]{\ensuremath{\mathbf{#1}}}
\newcommand*{\dfequ}[1][df]{\stackrel{\text{#1}}{=}}
\newcommand*{\bra}[1]{\ensuremath{\langle #1 \rvert}}
\newcommand*{\ket}[1]{\ensuremath{\lvert #1 \rangle}}
\newcommand*{\braket}[2]{\ensuremath{\langle #1 | #2 \rangle}}


%% =====================================================================
%%

\begin{document}

\newcounter{seria}
\setcounter{seria}{6}
% 1 - sortowanie
% 2 - układy równań
% 3 - układy równań: iteracyjnie
% 4 - wartości własne metodą potęgową
% 5 - wartości własne przez diagonalizację
% 6 - równania różniczkowe zwyczajne: Runge-Kutta
\newcommand*{\preparedby}{Jakub Narębski}

\title{
  \bfseries
  Numeryczne rozwiązywanie równań różniczkowych zwyczajnych cz. I\\
  \normalfont
  Metody typu Runge--Kutty
  \\[1ex]
  \Large
  Ćwiczenia z~Metod Numerycznych II~A\\
  rok akademicki 2001/2002\\[0.5ex]
  \normalsize seria \arabic{seria}}

\author{
  \href{mailto:jnareb@fuw.edu.pl}{Jakub Narębski} \and
  \href{mailto:tpawlow@fuw.edu.pl}{Tomasz Pawłowski}}

\date{\today}

\maketitle


\section{Wstęp}
\label{sec:wstep}

Rozwiązywanie równań różniczkowych zwyczajnych można zawsze sprowadzić
do rozwiązywania układu równań różniczkowych pierwszego rzędu.
Najprościej jest wybrać nowe zmienne tak, by były po prostu pochodnymi
tzn. $\vect{u} = [u(t), u'(t), u''(t), \dots]$.  Wystarczy zatem
rozważać układ równań różniczkowych pierwszego rzędu.  Rozpatrywać
będziemy \emph{zagadnienia początkowe}, tzn. będziemy chcieli znaleźć
rozwiązanie równania różniczkowego dla~zadanej wartości \(\vect{u}_0\)
w~punkcie początkowym \(t_0\):
\begin{equation}
  \label{eq:ODE}
  \left\{
  \begin{aligned}
    %\label{eq:ODE:eq}
    \frac{d \vect{u}(t)}{d t} &= \vect{f}(\vect{u},t)\\
    %\label{eq:ODE:start}
    \vect{u}(t_0) &= \vect{u}_0
  \end{aligned}
  \right.
\end{equation}

Rozważać będziemy \emph{metody różnicowe}, które wyznaczają
przybliżenia jedynie wartości rozwiązania \(\vect{u}(t_i)\) w~pewnych
punktach \(t_i\).  Będziemy się zajmowali przede wszystkim przypadkiem,
gdy punkty \(t_i\) są rozmieszczone równomiernie co \(\Delta t\),
tzn. \(t_i = t_0 + i\, \Delta t\).  Wielkość \(\Delta t\) nazywamy
\emph{krokiem całkowania}.  Metody różnicowe dają nam wzór na
wyliczenie wartości funkcji w~kroku następnym \(t_{n+1}\) na podstawie
kroków~poprzednich.

Dla uproszczenia wprowadzimy następujące oznaczenia
\begin{subequations}
  \label{eq:oznaczenia}
  \begin{align}
    \label{eq:oznaczenia:t_i}
    t_i &= t_0 + i\: \Delta t \\
    \label{eq:oznaczenia:u_i}
    \vect{u}_i &= \vect{u}(t_i) \\
    \label{eq:oznaczenia:f_i}
    \vect{f}_i &= \vect{f}(\vect{u}_i, t_i)
  \end{align}
\end{subequations}

Ogólną metodą na wyprowadzenie metod różnicowych rozwiązywania układu
równań różniczkowych zwyczajnych jest zapisanie go w~postaci całkowej
(poprzez odcałkowanie równania różniczkowego po jednym kroku czasowym,
od \(t_n\) do \(t_{n+1}\)) i~zastosowanie przybliżenia na całkę po
prawej stronie (całkę z~funkcji \(\vect{f}(\vect{u}(t),t)\)):
\begin{subequations}
  \begin{align*}
    \frac{d\vect{u}}{dt} &= \vect{f}(\vect{u}(t),t) \\
    \int_{t_n}^{t_{n+1}} \frac{d\vect{u}}{dt}\: dt
    &= \int_{t_n}^{t_{n+1}} \vect{f}(\vect{u}(t),t)\: dt \\
    \vect{u}(t_{n+1}) - \vect{u}(t_n)
    &= \int_{t_n}^{t_{n+1}} \vect{f}(\vect{u}(t),t)\: dt 
  \end{align*}
\end{subequations}


\section{Metoda Eulera}
\label{sec:metoda-eulera}

Najprostszą metodą różnicową jest \emph{metoda Eulera}.  Pochodną
\(\vect{u}'(t)\) przybliżamy ilorazem różnicowym pierwszego rzędu
w~przód opartym na węzłach \(t_n\) i~\(t_{n+1} = t_n + \Delta
t\).  Inaczej: przybliżamy całkę za pomocą wartości funkcji w~początku
przedziału razy jego długość:
\begin{equation*}
  \int_{t_n}^{t_{n+1}} \vect{f}(\vect{u}(t),t)\: dt 
  \simeq \Delta t \, \vect{u}(t_n)
\end{equation*}

Dostajemy zatem następujący wzór na wartość rozwiązania w~kroku
następnym
\begin{equation}
  \label{eq:Euler}
  \vect{u}_{n+1} = \vect{u}_n + \Delta t\, \vect{f}(\vect{u}_n,t_n)  
\end{equation}

Metoda ta ma prostą interpretację geometryczną.  Rozpoczynając od
punktu \((t_i,\vect{u}_i)\) leżącego na krzywej całkowej będącej
rozwiązaniem zagadnienia~\eqref{eq:ODE} posuwamy się po stycznej do
wykresu aż osiągniemy punkt o~współrzędnej niezależnej \(t_{i+1}\).




\section{Ulepszenia metody Eulera}
\label{sec:ulepsz-metody-Eulera}

Metoda Eulera jest rzadko stosowana ze względu na swoją bardzo wolną
zbieżność (musimy stosować bardzo mały krok by mieć dobrą dokładność)
i~niestabilność (w~niektórych zagadnieniach dla dowolnie małego kroku
nie uzyskamy dobrego przybliżenia).  Można ją ulepszyć, biorąc styczną
do krzywej nie w~punkcie \(t_i\) tylko w~połowie przedziału,
t.j. w~\(t_{i+1/2}\).  Wartość funkcji w~połowie przedziału~wyliczamy
metodą Eulera z~dwukrotnie mniejszym krokiem:
\begin{subequations}
  \label{eq:ulepszony-Euler}
  \begin{align}
    \label{eq:ulepszony-Euler:u-half}
    \vect{u}_{n+1/2} 
    &= \vect{u}_n + \frac{\Delta t}{2}\, \vect{f}(\vect{u}_n,t_n) \\
    \label{eq:ulepszony-Euler:u-next}
    \vect{u}_{n+1\hphantom{/2}} 
    &= \vect{u}_n + \Delta t\, \vect{f}(\vect{u}_{n+1/2},t_{n+1/2}) 
  \end{align}
\end{subequations}
co można zapisać w~skróconej postaci
\begin{equation}
  \label{eq:ulepszony-Euler:skompresow}
  %\tag{\ref{eq:ulepszony-Euler}'}
  \vect{u}_{n+1} = \vect{u}_n 
  + \Delta t\:\vect{f}(\vect{u} + \frac{\Delta t}{2}\,\vect{f}(\vect{u}_n,t_n), t_{n+1/2}) 
\end{equation}
Tą metodę (te ulepszenie metody Eulera) nazywa się \emph{metodą midpoint}
lub \emph{zmodyfikowaną metodą Eulera}.

Można także wziąć inne ulepszenie, tzn. jako przybliżenie ilorazu
różnicowego biorąc średnią arytmetyczną stycznych w~punktach \(t_n\)
i~\(t_{n+1}\).  Wartość stycznej w~punkcie \(t_{n+1}\) obliczamy za
pomocą zwykłej metody Eulera:
\begin{subequations}
  \label{eq:Heun}
  \begin{align}
    \label{eq:Heun:k1}
    \vect{k}
    &= \vect{u}_n + \Delta t\,\vect{f}(\vect{u}_n, t_n) \\
    \label{eq:Heun:u-next}
    \vect{u}_{n+1} 
    &= \vect{u}_n + \frac{1}{2} \Delta t\,
    \bigl(\vect{f}(\vect{u}_n, t_n) + \vect{f}(\vect{k}, t_{n+1})\bigr)
  \end{align}
\end{subequations}
co można zapisać w~skróconej postaci
\begin{equation}
  \label{eq:Heun:skompresow}
  %\tag{\ref{eq:Heun}'}
  \vect{u}_{n+1} 
  = \vect{u}_n + \frac{1}{2} \Delta t\,
  \bigl(
    \vect{f}(\vect{u})n, t_n) + 
    \vect{f}(\vect{u}_n + \Delta t\,\vect{f}(\vect{u}_n, t_n), t_{n+1})
  \bigr)
\end{equation}
Metodę tą nazywamy ,,ulepszoną metodą Eulera'' lub \emph{metodą Heuna}.



\section{Metoda Runge--Kutty}
\label{sec:metoda-runge-kutty}

Dzięki użyciu wartości funkcji pośrodku przedziału (ogólnie: dzięki
uśrednieniu) w~ulepszeniach metody Eulera osiągamy szybszą zbieżność
metody.  Skoro uśrednianie może pomóc, spróbujmy zastosować kilka
sposobów wyliczenia wartości funkcji pośrodku przedziału i~uśrednijmy
je (zastosujmy je wszystkie z~odpowiednimi wagami).
\begin{subequations}
  \label{eq:to-RK}
  \begin{align}
    \label{eq:to-RK:1}
    \frac{\vect{u}_{n+1} - \vect{u}_n}{\Delta t} &=
    \frac{1}{2} \bigl[
      \vect{f}(\vect{u}_n,t) + \vect{f}(\vect{u}_{n+1},t)
    \bigr] \\
    \label{eq:to-RK:2}
    \frac{\vect{u}_{n+1} - \vect{u}_n}{\Delta t} &=
    \vect{f} \Bigl(
      \vect{u}_n + \frac{\Delta t}{2}\vect{f}(\vect{u}_n,t_n),
      t_n + \frac{\Delta t}{2}
    \Bigr) \\
    \label{eq:to-RK:3}
    \frac{\vect{u}_{n+1} - \vect{u}_n}{\Delta t} &=
    \vect{f} \Bigl(
      \frac{\vect{u}_n + \vect{u}_{n+1}}{2},
      t_n + \frac{\Delta t}{2}
    \Bigr)
  \end{align}
\end{subequations}

Wprowadźmy następujące oznaczenia (uwaga: \(k_i\) są tutaj wektorami):
\begin{subequations}
  \label{eq:pre-RK:kn}
  \begin{align}
    \label{eq:pre-RK:k1}
    k_1 &= \vect{f} (\vect{u}_n, t_n) \\
    \label{eq:pre-RK:k2}
    k_2 &= \vect{f} (\vect{u}_{n+1/2}, t_{n+1/2}) \\
    \label{eq:pre-RK:k3}
    k_3 &= \vect{f} 
    \Bigl(\frac{\vect{u}_n + \vect{u}_{n+1}}{2},
    t_{n+1/2}\Bigr) \\
    \label{eq:pre-RK:k4}
    k_4 &= \vect{f}(\vect{u}_{n+1},t_{n+1})
  \end{align}
\end{subequations}

Aby móc efektywnie wyliczać \(k_i\) musimy wyrazić je przez wielkości
już wyliczone, t.j. wyrazić \(k_n\) przez \(k_i, i=1\dots n-1\).
Przekształcamy równanie~\eqref{eq:pre-RK:k3} by pozbyć się
\(\vect{u}_{n+1}\) po prawej stronie
\begin{equation*}
  \begin{split}
    \vect{f} 
    \Bigl(
    \frac{\vect{u}_n + \vect{u}_{n+1}}{2},
    t_n + \frac{\Delta t}{2} 
    \Bigr) 
    &=
    \vect{f} 
    \Bigl(
    \frac{\vect{u}_n + \vect{u}_n + (\vect{u}_{n+1} - \vect{u}_n)}{2},
    t_n + \frac{\Delta t}{2}
    \Bigr) 
    =\\
    &=
    \vect{f}
    \Bigl(
    \vect{u}_n + 
    \frac{\Delta t}{2} \frac{\vect{u}_{n+1} - \vect{u}_n}{\Delta t},
    t_n + \frac{\Delta t}{2}
    \Bigr)
    =\\
    &=
    \vect{f}
    \Bigl(
    \vect{u}_n + \frac{\Delta t}{2} k_2,
    t_n + \frac{\Delta t}{2}
    \Bigr)
  \end{split}
\end{equation*}
Skorzystaliśmy tutaj z~przybliżenia \(k_2\) na pochodną (iloraz
różnicowy).

Analogicznie pozbywamy się zależności od \(\vect{u}_{n+1}\) po prawej
stronie w~definicji \(k_4\), tym razem korzystając z~przybliżenia
\(k_3\) na pochodną (a~właściwie na iloraz różnicowy).

Dostajemy zatem
\begin{subequations}
  \label{eq:RK:kn}
  \begin{align}
    \label{eq:RK:k1}
    k_1 &= \vect{f}      ( \vect{u}_n                         , t_n ) \\
    \label{eq:RK:k2}
    k_2 &= \vect{f} \Bigl( \vect{u}_n + \frac{\Delta t}{2} k_1, t_n + \frac{\Delta t}{2} \Bigr) \\ 
    \label{eq:RK:k3}
    k_3 &= \vect{f} \Bigl( \vect{u}_n + \frac{\Delta t}{2} k_2, t_n + \frac{\Delta t}{2} \Bigr) \\ 
    \label{eq:RK:k4}
    k_4 &= \vect{f}      ( \vect{u}_n + \Delta t\, k_3        , t_n + \Delta t) 
  \end{align}
\end{subequations}

Aby wyliczyć dobre przybliżenie na pochodną (na iloraz różnicowy)
uśredniamy wzory~(\ref{eq:to-RK:1}\dots\ref{eq:to-RK:3}) (biorąc je
z~równymi wagami).
\begin{equation}
  \label{eq:pre-RK}
  \frac{\vect{u}_{n+1} - \vect{u}_n}{\Delta t} =
  \frac{1}{3} \bigl( \eqref{eq:to-RK:1} + \eqref{eq:to-RK:2} + \eqref{eq:to-RK:3} \bigr) =
  \frac{1}{6} \bigl( k_1 + 2 k_2 + 2 k_3 + k_4 \bigr)
\end{equation}

Ostatecznie dostajemy następującą metodę (klasyczną metodę
Runge--Kutty czwartego rzędu)
\begin{equation}
  \label{eq:RK}
  \vect{u}_{n+1} = \vect{u}_n +   \frac{\Delta t}{6} ( k_1 + 2 k_2 + 2 k_3 + k_4 )
\end{equation}

Warto zwrócić uwagę, że w~metodzie Runge--Kutty musimy wielokrotnie
obliczać wartość funkcji.  Jeśli zatem obliczanie wartości funkcji
zajmuje dużo czasu, to lepiej stosować inne metody.  Z~kolei
w~metodzie tej korzystamy tylko z~poprzedniego kroku, tak więc w~łatwy
sposób można za pomocą ekstrapolacji Richardsona zrobić z~niej metodę
o~zmiennym kroku.

Istnieją także inne wersje metody Rungego--Kutty czwartego rzędu,
gdzie bierze się wartości w~innych punktach i~składa je z~innymi
wagami.  Wzór klasyczny wykorzystuje się najczęściej ze względu na
jego prostotę.  Wzór Ralstona daje najlepsze oszacowanie błędu.  Przy
stosowaniu wzoru Gilla obliczenia można dokonywać tak, by wyrównywać
powstające błędy zaokrągleń.  Wzór ,,trzech ósmych'' jest niemal
równie prosty jak klasyczny; każdy z~\(k_i\) wyznacza wartość w~innym
punkcie: \(t_n, t_{n+1/3}, t_{n+2/3}, t_{n+1}\).  Istnieje także wzór
zamknięty na metodę Runge--Kutty piątego rzędu, który wykorzystuje się
przy konstrukcji metody ze zmiennym krokiem.

\pagebreak[4]

\section{Podsumowanie}
\label{sec:podsumowanie}

\parindent=0pt

Zmodyfikowana metoda Eulera
\begin{subequations}
  \label{eq:RK2mod-Euler}
  \begin{align}
    \label{eq:RK2mod-Euler:k1}
    k_1 &= \vect{f} (\vect{u}_n                           , t_n) \\
    \label{eq:RK2mod-Euler:k2}
    k_2 &= \vect{f} (\vect{u}_n + \frac{\Delta t}{2}\, k_1, t_n + \frac{\Delta t}{2}) \\[0.5ex]
    \label{eq:RK2mod-Euler:u_next}
    \vect{u}_{n+1} &= \vect{u}_n + \Delta t\, k_2 
  \end{align}
\end{subequations}

Metoda Heuna (ulepszona metoda Eulera)
\begin{subequations}
  \label{eq:RK2Heun}
  \begin{align}
    \label{eq:RK2Heun:k1}
    k_1 &= \vect{f} (\vect{u}_n                 , t_n) \\
    \label{eq:RK2Heun:k2}
    k_2 &= \vect{f} (\vect{u}_n + \Delta t\, k_1, t_n + \Delta t)\\[0.5 ex]
    \label{eq:RK2Heun:u_next}
    \vect{u}_{n+1} &= \vect{u}_n + \frac{\Delta t}{2} ( k_1 + k_2 )
  \end{align}
\end{subequations}

Wzór klasyczny
\begin{subequations}
  \label{eq:RK4class}
  \begin{align}
    \label{eq:RK4class:k1}
    k_1 &= \vect{f}      ( \vect{u}_n                         , t_n ) \\
    \label{eq:RK4class:k2}
    k_2 &= \vect{f} \Bigl( \vect{u}_n + \frac{\Delta t}{2} k_1, t_n + \frac{\Delta t}{2} \Bigr) \\ 
    \label{eq:RK4class:k3}
    k_3 &= \vect{f} \Bigl( \vect{u}_n + \frac{\Delta t}{2} k_2, t_n + \frac{\Delta t}{2} \Bigr) \\ 
    \label{eq:RK4class:k4}
    k_4 &= \vect{f}      ( \vect{u}_n + \Delta t\, k_3        , t_n + \Delta t) \\[0.5 ex]
    %%
    \label{eq:RK4class:u_next}
    \vect{u}_{n+1} &= \vect{u}_n +   \frac{\Delta t}{6} ( k_1 + 2 k_2 + 2 k_3 + k_4 )
  \end{align}
\end{subequations}


\noindent
Wzór ,,trzech ósmych''
\begin{subequations}
  \label{eq:RK4-3/8}
  \begin{align}
    \label{eq:RK4-3/8:k1}
    k_1 &= \vect{f}      ( \vect{u}_n                         , t_n ) \\
    \label{eq:RK4-3/8:k2}
    k_2 &= \vect{f} \Bigl( \vect{u}_n + \frac{1}{3}\Delta t\, k_1, 
                           t_n + \frac{1}{2}\Delta t \Bigr) \\ 
    \label{eq:RK4-3/8:k3}
    k_3 &= \vect{f} \Bigl( \vect{u}_n - \frac{1}{3}\Delta t\, k_1 + \Delta t\, k_2, 
                           t_n + \frac{2}{3}\Delta t \Bigr) \\ 
    \label{eq:RK4-3/8:k4}
    k_4 &= \vect{f}      ( \vect{u}_n + \Delta t\, k_1 - \Delta t\, k_2 + \Delta t\, k_3, 
                           t_n + \Delta t) \\[0.5 ex]
    %%
    \label{eq:RK4-3/8:u_next}
    \vect{u}_{n+1} &= \vect{u}_n +   \frac{\Delta t}{8} ( k_1 + 3 k_2 + 3 k_3 + k_4 )
  \end{align}
\end{subequations}


%
% USUNIĘTE
%
\begin{comment}

\section{Implementacja klasycznej metody Runge--Kutty czwartego rzędu}
\label{sec:implement-RK4}

Jedna z~implementacji klasycznej metody Runge--Kutty czwartego
rzędu.  Rozwiązuje układ równań różniczkowych zwyczajnych zadanych za
pomocą wektora funkcyj prawych stron.

\lstset{
  language=C,
  tabsize=2,
  indent=0pt,
  extendedchars=true
}

\small
\begin{lstlisting}{RK4-JN}
typedef double real;
#define REAL "%lg"
#define _R "l"

#define ZERO 2.0E-150

#define sqr(x) ((x)*(x))

#define FALSE 0
#define TRUE 1

#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))

/* czesc prawej strona rownania rozniczkowego zwyczajnego */
/* pierwszego rzedu                                       */
/*                                                        */
/*  x1'  = f1( t, x1,x2...,xN )                           */
/*  x2'  = f2( t, x1,x2...,xN )                           */
/*  ...                                                   */
/*  xN'  = fN( t, x1,x2...,xN )                           */
/*                                                        */
typedef real (*p_dequation1)( real t, real* x);

void RungeKutta4_classic( 
  real *res,       /* wynik: u_n+1  */
  int dim,         /* wymiar ukladu */
  p_dequation1 *f, /* prawa strona  */
  real* x,         /* punkt u_n     */
  real t,          
  real h )         /* krok: Dt      */
{
  int i;
  real *k1,*k2,*k3,*k4;
  real *tmp;

  /* UWAGA: powinno sie sprawdzac, czy alokacja pamieci */
  /*        sie udala czy nie                           */
  k1  = (real *)calloc(dim,sizeof(real));
  k2  = (real *)calloc(dim,sizeof(real));
  k3  = (real *)calloc(dim,sizeof(real));
  k4  = (real *)calloc(dim,sizeof(real));
  tmp = (real *)calloc(dim,sizeof(real));

  for(i = 0;i < dim; i++)
    k1[i] = (*f[i])(t,x);
  for(i = 0;i < dim; i++)
    tmp[i] = x[i] + 0.5*h*k1[i];
  for(i = 0;i < dim; i++)
    k2[i] = (*f[i])(t + 0.5*h,tmp);
  for(i = 0;i < dim; i++)
    tmp[i] = x[i] + 0.5*h*k2[i];
  for(i = 0;i < dim; i++)
    k3[i] = (*f[i])(t + 0.5*h,tmp);
  for(i = 0;i < dim; i++)
    tmp[i] =  x[i] + h*k3[i];
  for(i = 0;i < dim; i++)
    k4[i] = (*f[i])(t + h    ,tmp);

  for(i = 0;i < dim; i++)
    res[i] = x[i] + h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;

  free(tmp);
  free(k4);
  free(k3);
  free(k2);
  free(k1);
}
\end{lstlisting}

\normalsize
\bigskip
Implementacja metody Runge--Kutty z~,,Numerical Receipies in C''

\small
\begin{lstlisting}{RK4-NRC}
#define NRANSI
#include "nrutil.h"

void rk4(float y[], float dydx[], int n, float x, float h, float yout[],
         void (*derivs)(float, float [], float []))
/* Majac dane wartosci zmiennych y[1..n] i ich pochodnych dxdy[1..n]   */
/* w punkcie x, stosuje metode Runge-Kutty do obliczenia rozwiazania   */
/* w punkcie x+h. Zwraca obliczone rozwiazanie w yout[1..n], ktory nie */
/* musi byc rozlaczny z y (t.j. moze byc to ta sama tablica).          */
/* Uzytkownik powinien dostarczyc funkcje derivs(x,y,dxdy), ktora      */
/* zwraca pochodne dxdy w punkcie x                                    */
{
  int i;
  float xh,hh,h6,*dym,*dyt,*yt;

  /* alokacja pamieci */
  dym=vector(1,n);
  dyt=vector(1,n);
  yt=vector(1,n);

  hh=h*0.5;
  h6=h/6.0;
  xh=x+hh;

  /* wyliczenie k_1,...,k_4 */
  for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i];
    (*derivs)(xh,yt,dyt);
  for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i];
    (*derivs)(xh,yt,dym);
  for (i=1;i<=n;i++) {
    yt[i]=y[i]+h*dym[i];
    dym[i] += dyt[i];
  }
  (*derivs)(x+h,yt,dyt);

  /* wyliczenie rozwiazania */
  for (i=1;i<=n;i++)
    yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);

  /* zwolnienie pamieci */
  free_vector(yt,1,n);
  free_vector(dyt,1,n);
  free_vector(dym,1,n);
}
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software ]*!. */
\end{lstlisting}

\end{comment}

%%
\input{podpis.tex}

\end{document}



