[Pascal] Równanie różniczkowe metodą numeryczną

[galaj666]

mam problem z rozwiązaniem rownania rożniczkowego metodą numeryczna w jezyku pascal.
zadanie jakie otrzymałem to:
W reaktorze przepływowym o stałej temperaturze koncentrację składnika określa równianie
\(\displaystyle{ Q \frac{dc}{dV}=-R(c), c(0)=c_{in},
R(c)= \frac{k c}{c+v}}\)
gdzie \(\displaystyle{ Q=50, c_{in}=2, k=0,198, v=3,8}\). znaleźć zależnośc stęzenia składnika od bezwymiarowej objętości \(\displaystyle{ V}\) w przedziale \(\displaystyle{ [0, 2000]}\)

[SlotaWoj]

Przy numerycznym rozwiązywaniu równań różniczkowych przy pomocy komputera od samych równań ważniejsza jest używana (zaimplementowana) metoda, a ta może być różnego rzędu (im wyższy rząd tym większa dokładność i stabilność), jawna lub niejawna, stałokrokowa lub zmiennokrokowa (adaptacyjna), bez korekcji, predykcyjna, korekcyjna lub kombinowana.

Czy wydający temat zadania określił jaka metoda rozwiązania ma być użyta?
W grę wchodzą (jako najprostsze) metody: Eulera (jawna i niejawna) oraz Rungego-Kutty 2-go rzędu.

Generalnie numeryczne rozwiązywanie równań różniczkowych odbywa się iteracyjnie, tzn. na podstawie pochodnej i wartości rozwiązania uzyskanej w kroku \(\displaystyle{ n}\) (poprzednim) oblicza się wartość w kroku \(\displaystyle{ n+1}\) (bieżącym).

Przekształca się układ równań różniczkowych do postaci:

• \(\displaystyle{ \vec{y}(x)'=f(\vec{y}(x);x)}\)

gdzie \(\displaystyle{ y(x)}\) to rozwiązanie (funkcja), które trzeba wyznaczyć, a \(\displaystyle{ x}\) to zmienna niezależna (np. czas).

Gdy jest jedno równanie różniczkowe pierwszego rzędu, to ww. postać jest skalarna.

Dla najprostszej metody (jawnej Eulera) iteracje są zdefiniowane w sposób następujący:

• \(\displaystyle{ \vec{y}(x_{n+1})=f(\vec{y}(x_n);x_n)\mbox{·}h+\vec{y}(x_n)}\)

\(\displaystyle{ x_0}\) i \(\displaystyle{ \vec{y}(x_0)}\) są warunkami początkowymi, a \(\displaystyle{ h=x_{n+1}-x_n}\) krokiem iteracji.

W Twoim przypadku będzie:

• \(\displaystyle{ \frac{dc(V)}{dV}=f(c(V);V)=\frac{-kc(V)}{Q\left(c(V)-v\right)} \\
  V_0=0 \\
  c(V_0)=c_{in}}\)

Zwróć uwagę, że \(\displaystyle{ \frac{-kc(V)}{Q\left(c(V)-v\right)}}\) nie zależy bezpośrednio od \(\displaystyle{ V}\).

Gdy przyjmiemy \(\displaystyle{ h=\Delta V=20}\) to w przedziale \(\displaystyle{ \left[0;2000\right]}\) będzie 100 iteracji.

Gdy uzyskana w ten sposób wartość \(\displaystyle{ c(2000)}\) będzie różnić się znacząco od analogicznej wartości uzyskanej np. dla \(\displaystyle{ h=10}\) (200 iteracji) oznacza, że krok \(\displaystyle{ h=20}\) był za duży i wynik uzyskany dla kroku \(\displaystyle{ h=10}\) należy porównać z wynikiem dla kroku jeszcze mniejszego (np. \(\displaystyle{ h=5}\)).-- 16 kwi 2015, o 16:20 --Zarys kodu programu z zaimplementowana jawną metodą Eulera będzie wyglądał następująco:

Kod: Zaznacz cały

Program  RRREuler;
  { Deklaracje }
Function F (c, V : extended ) : extended;  { extended: na PC najszerszy typ float }
  begin
     F := -k*c/Q/(c-vvv)  { vvv=3,8 }
  end;
begin
  { Inicjacja i wczytanie danych }
  h := cmax/itenum;    { itenum: liczba kroków }
  V := V0;  c := cin;  { V0=0 }
  for  i = 1 to itenum
    do
      begin
         nextV := F(c,V)*h+V;  { tu można bezpośrednio przypisać do V }
         c := c + h;
         V := nextV;           { patrz komentarz powyżej }
      end;
  { Zakończenie }
end.

Oczywiście deklaracje trzeba rozpisać, a program rozbudować o to co „zakomentowałem”, oraz np. o ocenę dokładności etc.

Oczywiście jądro programu można zakodować jako procedurę, aby można było ją wykonywać wielokrotnie, a podczas iteracji wypełniać tablicę z kolejnymi wartościami rozwiązania, aby można było wykonać wykres, etc.

Jeśli znasz choć trochę Pascal i potrafisz ułożyć proste algorytmy przetwarzające np. tablice jednowymiarowe (szukanie minimum, liczenie elementów należących do pewnego przedziału, sortowanie, etc.) to krok po kroku osiągniesz zadowalający efekt.