Numeryczne rozwiązywanie równań różniczkowych zwyczajnych.

Równania różniczkowe opisują wiele problemów technicznych. Przebieg zjawisk można często doprowadzić do postaci układu równań pierwszego rzędu o postaci:

, (1)

z danymi warunkami początkowymi u(t₀)=u₀.

Dla układu N równań u, u₀ są wektorami o N elementach, a f zbiorem funkcji.

Często rozwiązania tego układu nie możemy znaleźć analitycznie. Musimy wtedy zadowolić się rozwiązaniem numerycznym. Rozwiązanie to stanowi zbiór wartości funkcji u(t) dla skończonej ilości wycinków należących do interesującego nas przedziału.

Najpowszechniej stosowaną metodą rozwiązywania równań różniczkowych dla problemów technicznych jest metoda Rungego-Kutty.

Metoda rozwiązania:

Zapiszmy rozwinięcie równania wyjściowego (1) w szereg Taylora w otoczeniu punktu u(t+h):

(2).

Utwórzmy teraz wyrażenia:

K₁=hf(t,u)

K₂=hf(t+b₂+h,u+l₂₁K₁)

K₃=hf(t+b₃+h,u+l₃₁K₁+ l₃₂K₂)

...........

Staramy się wyznaczyć stałe l_ij ,b_i, a_i, aby wyrażenie

, (3),

Było zgodne z prawą stroną rozwinięcia w szereg (2) z jak największą dokładnością. Korzystając z r elementów sumy (3) możemy uzyskać zgodność z wyrażeniami szeregu (2) do stopnia r względem h. Liczbowe wyznaczenie współczynników l, a,b dla metod wyższego rzędu jest skomplikowane. Dla czwartego rzędu współczynniki mogą mieć wartości:

a₁=a₄=1/6 a₂=a₃=1/3

b₁=0 b₂=b₃=1/2 b₄=1

l₂₁=1/2 l₃₁=0 l₃₂=1/2 l₃₁=0

l₄₁=l₄₂=0 l₄₃=1.

Równania wówczas przyjmują postać:

u_k+1=u_k+1/6(K₁+2K₂+2K₃+K₄)

K₁=hf(t_k,u_k)

K₂=hf(t_k+(1/2)h,u_k+(1/2)K₁)

K₃=hf(t_k+(1/2)h,u_k+(1/2)K₂)

K₃=hf(t_k+h,u_k+K₃).

Dla tej metody jednokrokowej należy obliczać r-krotnie wartość funkcji f w każdym kroku, tzn. wewnątrz przedziału (t_k,t_k+h).

Biblioteka ruku.tpu zawierająca procedurę:

procedure Runge_Kutta_4( liczba_rownan :integer;

t0 :real;

krok :real;

var u :wektor;

układ_rownan :procedura);

korzystając z wyżej pokazanych wzorów wyznacza wartości szukanych funkcji i zawiera je w wektorze u. Są to wartości przybliżone i w punkcie t_o+krok. Na wejściu parametr u zawiera wartości początkowe zmiennych niezależnych.

Procedura:

procedure układ_rowanan (t:real; u:wektor; var du:wektor);

t - wartość zmiennej niezależnej,

u - wartości wektora zmiennych zależnych,

du - obliczone prawe strony układu równań w punkcie (t,u) (wartości funkcji f_i).

Należy zwrócić uwagę, że procedura Runge_Kutte_4 oblicza jedynie jeden punkt rozwiązania. W celu wyznaczenia rozwiązania różniczkowego w jakimś przedziale należy wywołać tą procedurę cyklicznie wprowadzając wyniki poprzedniego kroku jako wartości początkowe.

Załączony program o nazwie pocisk_Robert_Szmurlo symuluje lot pocisku o zadanej masie, prędkości początkowej i kącie alfa. W celu zbliżenia toru do warunków rzeczywistych program bierze również pod uwagę współczynnik oporu powietrza.

Procedury:

skalowanie służy do wyskalowania obrazu na ekranie w zależności od podanych parametrów.

rysowanie_tla generuje linie osi oraz opisy do nich.

rysowanie_wykresu sparametryzowana procedura rysująca wykres toru w zależności od zadanych wartości.

wczytaj_wartosc procedura wczytywania wartości do której przekazuje się wartość początkową, wartość minimalną i maksymalną. Po wciśnięciu samego entera do zmiennej przekazywa jest wartość początkowa.

część_główna Część główna programu steruje całym cyklem programu. Po przeanalizowaniu wykresu użytkownik może dokonać zmian wybranych parametrów. Następnie wykres jest ponownie generowany z nowymi parametrami.

Sparametryzowanie procedury rysowania wykresu oraz skalowania umożliwia szybkie zwiększenie ilości wykresów, widocznych jednocześnie. W tym celu niezbędne poprawki wymagane są jedynie w głównej części programu bez konieczności zmian procedur.

Konieczność rozwiązywania równań różniczkowych nie ogranicza się do wyznaczenia toru pocisku. Z równaniami różniczkowymi stykają się wszyscy inżynierowie. Równania te opisują zasadniczą większość procesów technicznych. Rozwiązanie problemów elektrycznych, elektronicznych w stanach nieustalonych jest praktycznie niemożliwe bez rozwiązania równań różniczkowych.

PROGRAM pocisk_Robert_Szmurlo;

USES crt,ruku,graph;

CONST

g=9.81;

VAR

b,m,m1 :real;

tp,tk,h,tymax,wys,wys1 :real;

x,y :wektor;

zas,alfa,Vo,zas1,alfa1,Vo1 :real;

l1,l2,l :real;

k :char;

q,r,czas,li :string [31];

sterownik,tryb,i,wysmax :integer;

PROCEDURE uklad_row_y (t: real; u:wektor; var du:wektor); far;

begin

du [1] := u[2];

du [2] := - b/m * u[2] - g ;

end;

PROCEDURE uklad_row_x (t: real; u:wektor; var du:wektor); far;

begin

du [1] := u[2];

du [2] := - b/m * u[2] ;

end;

PROCEDURE skalowanie(vpocz,masa,kat:real;var krok:real);

var

krokx,kroky :real;

a :integer;

begin

a:=trunc((vpocz+masa)/2);

case a of

5001..10000 :h:=10;

3001..5000 :h:=5;

2001..3000 :h:=2;

1501..2000 :h:=1;

1001..1500 :h:=0.5;

501..1000 :h:=0.1;

0..500 :h:=0.01

else h:=1;

end;

tp:=0;

x[1]:=0;x[2]:=vpocz*cos(kat*pi/180);

y[1]:=0;y[2]:=vpocz*sin(kat*pi/180);

wys:=0;

while y[1]>=0 do

begin

Runge_Kutta_4 (2, tp, h,x, uklad_row_x);

Runge_Kutta_4 (2, tp, h,y, uklad_row_y);

if (y[1]>wys) then wys:=y[1];

tp:=tp+h

end;

if x[1]=0 then x[1]:=0.000000001;

if wys=0 then wys:=0.000000001;

krokx:=580/x[1];

kroky:=300/wys;

if (kroky*x[1]>580) then krok:=krokx

else krok:=kroky;

end;

PROCEDURE rysowanie_tla;

begin

clearviewport;

line(10,400,620,400);

line(615,395,620,400);

line(615,405,620,400);

settextstyle(smallfont,horizdir,4);

outtextxy(550,410,'odleglosc [m]');

outtextxy(20,10,'wysokosc [m]');

line(10,10,10,400);

line(6,14,10,10);

line(14,14,10,10);

end;

procedure rysowanie_wykresu(vpocz,kat,masa,wsp:real;color:word;wspx,wspy,n:integer);

begin

wys:=0;

tp:=0;

x[1]:=0;x[2]:=vpocz*cos(kat*pi/180);

y[1]:=0;y[2]:=vpocz*sin(kat*pi/180);

moveto(10,400);

setcolor(color); {}

repeat

begin

Runge_Kutta_4 (2,tp,h,x,uklad_row_x);

Runge_Kutta_4 (2,tp,h,y,uklad_row_y);

if (y[1]>wys) then wys:=y[1];

lineto(10+round(l*x[1]),400-round(l*y[1]));

tp:=tp+h

end;

until y[1]<0;

wysmax:=400-round(wys*l);

line(8,wysmax,12,wysmax);

str(wys:30:1,r);

Runge_Kutta_4 (2, tp-h, h,x, uklad_row_x);

zas:=x[1];

str(zas:30:1,q);

str(tp-h:30:1,czas);

str(n:30,li);

line(10+round(zas*l),398,10+round(zas*l),402);

settextstyle(smallfont,horizdir,4);

outtextxy(wspx,wspy,'Maksymalna wysokosc'+r+' m');

outtextxy(wspx,wspy+10,'Zasieg rzutu '+q+' m');

outtextxy(wspx,wspy+20,'Czas lotu '+czas+' s');

end;

procedure wczytaj_wartosc(var wart:real;min,max,wart_pocz:real);

var

wart_moja :real;

xtext,ytext :byte;

slowo :string;

code :integer;

begin

xtext:=wherex;ytext:=wherey;

repeat

begin

gotoxy(1,1+ytext);

write('Wartosc dotychczasowa: ',wart_pocz:5:3,' Wartosc nowa: ');

readln(slowo);

if slowo='' then wart_moja:=wart_pocz

else val(slowo,wart_moja,code);

end;

until (wart_moja<=max)and (wart_moja>=min);

wart:=wart_moja;

end;

procedure inicjuj_dane;

begin

Vo:=0;Vo1:=0;alfa:=0;alfa1:=0;b:=0;m:=0;m1:=0;

end;

BEGIN

inicjuj_dane;

repeat

begin

clrscr;

textcolor(yellow);

write('Podaj predkosc poczatkowa Vo <100-2000>');wczytaj_wartosc(Vo,100,2000,Vo);

write('Podaj kat wyrzutu alfa <0-90> ');wczytaj_wartosc(alfa,0,90,alfa);

write('Podaj mase m<0.001-1000kg> ');wczytaj_wartosc(m,0.001,1000,m);

textcolor(cyan);

write('Podaj predkosc poczatkowa Vo1 <100-2000>');wczytaj_wartosc(Vo1,100,2000,Vo1);

write('Podaj kat wyrzutu alfa1 <0-90> ');wczytaj_wartosc(alfa1,0,90,alfa1);

write('Podaj mase m1<0.001-1000kg> ');wczytaj_wartosc(m1,0.001,1000,m1);

textcolor(white);

write('Podaj wspolczynnik oporu b>0 ');wczytaj_wartosc(b,0,10,b);

sterownik:=detect;

initgraph(sterownik,tryb,'c:\tp\bgi');

clearviewport;

outtextxy(230,230,'Jeszcze chwileczke ...');

skalowanie(Vo,m,alfa,l1);

skalowanie(Vo1,m1,alfa1,l2);

if l1>l2 then l:=l2

else l:=l1;

rysowanie_tla;

rysowanie_wykresu(Vo,alfa,m,b,yellow,120,10,1);

rysowanie_wykresu(Vo1,alfa1,m1,b,cyan,120,40,2);

setcolor(white);

settextstyle(triplexfont,horizdir,2);

outtextxy(150,410,'Rzut ukosny pocisku');

outtextxy(150,430,'Wyjscie z programu - "k"');

k:=readkey;

closegraph;

end;

until k='k';

END.

---