Informatyka I  Lab 09, r.a. 2011/2012
prow. SÅ‚awomir Czarnecki
Zadania na laboratorium nr. 9
Po utworzeniu nowego projektu, dołącz bibliotekę bibs.h.
1. Zadanie problemowe z równań ró\
niczkowych zwyczajnych  wstęp.
Algorytm Runge-Kutta rozwiązywania układu dwóch równań ró\
niczkowych zwyczajnych
pierwszego rzędu:
u ' x = f îÅ‚x,u x ,v x Å‚Å‚ ,
( ) ( ) ( )ûÅ‚
ðÅ‚
v ' x = g îÅ‚x,u x ,v x Å‚Å‚
( ) ( ) ( )ûÅ‚
ðÅ‚
z dwoma znanymi funkcjami
f : R3 R ; f = f x,u,v ,
( )
g : R3 R ; g = g x,u,v
( )
oraz z dwoma nieznanymi i poszukiwanymi funkcjami
u = u x , v = v x
( ) ( )
które spełniać mają warunki początkowe
u x0 = u0 , v x0 = v0 .
( ) ( )
Dane wejściowe i parametry algorytmu
" punkt poczÄ…tkowy (chwila poczÄ…tkowa x = x0): x0
" warunek poczÄ…tkowy dla funkcji u = u(x): u0
" warunek poczÄ…tkowy dla funkcji v = v(x): v0
" dÅ‚ugość kroku (czasowego): ¸
" liczba dyskretnych punktów (liczba chwil czasowych) n
w których obliczane bÄ™dÄ… wartoÅ›ci funkcji: u = u(i ·¸) , v = v(i ·¸ ) (i = 0,1,...,n  1)
" wskaz
nik do funkcji typu double f (double , double , double): *f
" wskaz
nik do funkcji typu double g (double , double , double): *g
" wektor punktów, w których obliczamy numerycznie rozwiązanie: X[n]
" wektor wartości funkcji u = u(x): U[n]
" wektor wartości funkcji v = v(x): V[n]
Algorytm Runge-Kutta  pseudo-kod
1.zainicjalizuj pierwsze składowe:
X[0] = x0
U[0] = u0
V[0] = v0
2. dla ka\
dej chwili czasowej i 0 d" i < n -1 , obliczaj
( )
{
F = ¸ · f (X[i] , U[i] , V[i])
G = ¸ · g (X[i] , U[i] , V[i])
dU = ¸ · f (X[i] + 0.5 ¸ , U[i] + 0.5 F , V[i] + 0.5 G)
dV = ¸ · g (X[i] + 0.5 ¸ , U[i] + 0.5 F , V[i] + 0.5 G)
X[ i + 1 ] = X[i] + ¸
U[ i + 1 ] = U[i] + dU
V[ i + 1 ] = V[i] + dV
}
Na podstawie powy\
szego algorytmu, zdefiniuj funkcjÄ™
void runge_kutta( double x0 , double u0 , double v0 , double ¸ , int n,
double(*f )(double x , double u , double v),
double(*g)(double x , double u , double v),
double* X , double* U , double* V )
i znajdz
następnie rozwiązanie numeryczne u = u t , v = v t układu dwóch równań
( ) ( )
ró\
niczkowych zwyczajnych:
u ' = 3u - v
v ' = 4u - v
z warunkami poczÄ…tkowymi
u 0 = 0.2
( )
v 0 = 0.5
( )
w przedziale T = 0 , 2 wywoÅ‚ujÄ…c funkcjÄ™ runge_kutta(& ) dla kroku czasowego ¸ = 0.05 .
[ ]
Sprawdz
poprawność rozwiązania numerycznego, konfrontując otrzymane wartości dyskretne
funkcji u = u t , v = v t z wartościami dyskretnymi, znanego w tym przypadku,
( ) ( )
et t et et t et
rozwiÄ…zania analitycznego: u t = - , v t = - .
( ) ( )
5 10 2 5
T
SporzÄ…dz
wykres (np. w Excel) funkcji wektorowej u,v t = îÅ‚u t ,v t Å‚Å‚ " !2 , który
( )( ) ( ) ( )ûÅ‚
ðÅ‚
powinien wyglądać tak jak na zamieszczonym poni\
ej rysunku rys.1.
1
0.8
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3
T
Rys.1. Wykres funkcji wektorowej u,v t = îÅ‚u t ,v t Å‚Å‚ " !2
( )( ) ( ) ( )ûÅ‚
ðÅ‚
et t et et t et
dla u t = - , v t = -
( ) ( )
5 10 2 5
Algorytm poszukiwania rozwiązania równania ró\
niczkowego zwyczajnego
rzędu drugiego
Zadanie numerycznego rozwiązania równania ró\
niczkowego zwyczajnego drugiego rzędu
y" x = h îÅ‚x, y x , y ' x Å‚Å‚
( ) ( ) ( )ûÅ‚
ðÅ‚
ze znanÄ… funkcjÄ…
h : R3 R
oraz nieznanÄ… i poszukiwanÄ… funkcjÄ…
y : R R
spełniającą warunki początkowe
y x0 = y0 , y ' x0 = v0
( ) ( )
gdzie x0 , y0 , v0 są znanymi wartościami, mo\
na z reguły odpowiednio przedefiniować w
terminach zadania rozwiązywania układu dwóch równań ró\
niczkowych zwyczajnych (x0 jest
najczęściej początkową chwilą czasu, y0 jest najczęściej wartością współrzędnej uogólnionej
Lagrange a definiującej poło\
enie punktu materialnego lub układu materialnego o jednym
stopniu swobody w chwili początkowej x0, natomiast v0 jest wtedy prędkością uogólnioną w
chwili poczÄ…tkowej x0).
Algorytm
" zdefiniuj funkcjÄ™
f : R3 R , f x,u,v = v
( )
oraz funkcjÄ™
g : R3 R , g x,u,v = h x,u,v
( ) ( )
gdzie u i v są odpowiednio wartościami funkcji:
u : R R, u x = y x i v : R R, v x = y ' x
( ) ( ) ( ) ( )
" Znajdz
rozwiązania u = u x i v = v x układu dwóch równań ró\
niczkowych
( ) ( )
zwyczajnych pierwszego rzędu
u ' x = f îÅ‚x,u x ,v x Å‚Å‚
( ) ( ) ( )ûÅ‚
ðÅ‚
v ' x = g îÅ‚x,u x ,v x Å‚Å‚
( ) ( ) ( )ûÅ‚
ðÅ‚
przy warunkach poczÄ…tkowych
u x0 = y0 , v x0 = v0
( ) ( )
wykorzystując na przykład opisany powy\
ej algorytm Runge-Kutta.
Zadanie problemowe  wahadło matematyczne
Rozwa\
my równanie ró\
niczkowe ruchu wahadłowego punktu materialnego o masie m,
poruszającego się na końcu niewa\
kiej linki o długości L, w polu grawitacyjnym ziemi
(a = 9.81 [m/s2]). Zakładamy dodatkowo, \
e istnieje siła tłumienia ruchu, która jest
proporcjonalna do prędkości kątowej ruchu linki. Współczynnik tłumienia oznaczmy przez c.
Przemieszczenie punktu w ka\
dej chwili czasu t określa kąt y = y(t). Prędkością kątową ruchu
obrotowego linki będzie zatem y = y (t).
2
d y c dy a dy
+ + sin y = 0, y 0 = Õ0 , 0 = 0
( ) ( ) ( )
dt2 mL dt L dt
Równanie ró\
niczkowe zapiszemy w postaci
2
d y a c dy
= - sin y -
( )
dt2 L mL dt
i definiujemy funkcje
f : R3 R ; f t,u,v = v
( )
a c
g : R3 R ; g t,u,v = - sin u - v
( ) ( )
L mL
Ponadto załó\
my, \
e
a = 9.81[m/s2] , m = 1.0 [kg] , L = 10.0 [m] , c = 3.0 [N s]
chwila poczÄ…tkowa: t0 = 0
Ä„ Ä„
warunek poczÄ…tkowy dla funkcji u = u(t): u0 = , tj. Õ0 = y t0 =
( )
2 2
dy
warunek poczÄ…tkowy dla funkcji v = v(t): v0 = 0 , tj. t0 = 0
( )
dt
dÅ‚ugość kroku (czasowego): ¸ = 0.001
liczba chwil czasowych: n = 30000
Znajdz
rozwiązanie numeryczne y = y t wywołując funkcję Runge_Kutta(& ), której
( )
wykres jest pokazany na rys.2.
2
1.5
1
0.5
0
0 10 20 30 40
-0.5
-1
-1.5
Rys. 2. Wykres funkcji y = y t dla wahadła matematycznego.
( )
Zadanie problemowe - oscylator harmoniczny
Rozwa\
my równanie ró\
niczkowe ruchu punktu materialnego o masie m, z warunkami
początkowymi na przemieszczenie y = y(t) i prędkość y = y (t)
2
d y
m + ky = F t
( )
dt2
dy
y 0 = 0 , 0 = 0
( ) ( )
dt
k m
F = F(t)
y = y(t)
Umieszczony na końcu sprę\
yny o stałej k punkt ma masę m. Do punktu materialnego
przyło\
ona jest, zale\
na tylko od czasu Ä, siÅ‚a F = F(Ä). Pomijamy tarcie pomiÄ™dzy punktem a
podło\
em.
Nietrudno wykazać, \
e całka ogólna tego równania (jego rozwiązanie) dana jest wzorem
t
1
y = y t =
( )
+"sin(Ét -ÉÄ ) F (Ä )dÄ , t e" 0 .
mÉ
0
gdzie zdefiniowano nową wielkość, tzw. częstość kołową drgań
k
É =
m
Dla siły F = F(t) danej wzorem
Fconst
Å„Å‚
Ä Ä " 0,tconst
[ ]
ôÅ‚
tconst
F Ä =
( )
òÅ‚
ôÅ‚
Fconst Ä e" tconst
ół
F(<)
Fconst
<
tconst
całkę ogólną równania ró\
niczkowego mo\
emy przedstawić w postaci następującego wzoru:
Å„Å‚
ëÅ‚ sin Ét öÅ‚
1 ( )
- t " 0,tconst
[ ]
ôÅ‚ ÷Å‚
tconst ìÅ‚t É
íÅ‚ Å‚Å‚
y t = ySTAT ôÅ‚
( )
òÅ‚
ôÅ‚1+ 1 îÅ‚sin Ét -Étconst - sin Ét Å‚Å‚ t e" tconst
( ) ( )ûÅ‚
ôÅ‚
Étconst ðÅ‚
ół
gdzie
Fconst Fconst
ëÅ‚ öÅ‚
ySTAT = =
ìÅ‚ ÷Å‚
k mÉ2 Å‚Å‚
íÅ‚
jest rozwiązaniem statycznym równania ró\
niczkowego, dla stałego w czasie obcią\
enia
Fconst.
W celu sprawdzenia skuteczności algorytmu Runge-Kutta rozwiązywania układu dwóch
równań ró\
niczkowych zwyczajnych, dopasujmy najpierw oznaczenia i zdefiniujmy
poprawnie wszystkie funkcje, które są parametrami naszej procedury. Najpierw równanie
2
d y
m + ky = F t
( )
dt2
zapiszmy w postaci
2
F t
d y ( ) - ky
=
dt2 m
i zgodnie z podanymi wzorami oraz po uwzględnieniu faktu, \
e rolÄ™ zmiennej x, odgrywa
teraz parametr t, mamy następujące definicje funkcji f i g:
f : R3 R ; f t,u,v = v
( )
g : R3 R ; g t,u,v = h t,u,v
( ) ( )
gdzie dla zdefiniowanej wcześniej funkcji F = F(t)
F t - ku
( )
h t,u,v =
( )
m
Wartości pozostałych parametrów są następujące:
chwila poczÄ…tkowa: x = 0
0
warunek poczÄ…tkowy dla funkcji u = u(t): u = 0
0
warunek poczÄ…tkowy dla funkcji v = v(t): v = 0
0
dÅ‚ugość kroku (czasowego): ¸ = 0.001
liczba chwil czasowych n, n = 500 ÷ 2000
Dla przyjętych danych, przeprowadz
symulację numerycznego rozwiązania naszego równania
oscylatora, dla kilku ró\
nych wartości chwili czasowej tconst charakteryzującej moment, w
którym narastająca liniowo siła F = F(t) osiąga stałą wartość Fconst = 1000.0 [N]. Przyjmij
ponadto stałą sprę\
yny k = m ·É2 dla m = 1.0 [kg] i É = 50.0 [1/s]. ChwilÄ™ czasowÄ… tconst
zdefiniuj następująco:
tconst = ¾ Å"T
gdzie
2Ä„
T =
É
natomiast ¾ jest dowolnÄ…, ale ustalonÄ… liczbÄ… dodatniÄ…. Wykresy sporzÄ…dz
nie dla funkcji y =
y(t), ale dla ilorazu
y t
( )
ySTAT
gdzie
1000.0 N
Fconst Fconst [ ]
ëÅ‚ öÅ‚
ySTAT = = = = 0.4[m] .
ìÅ‚
1
k mÉ2 ÷Å‚
íÅ‚ Å‚Å‚
1.0 kg Å"50.02 îÅ‚ Å‚Å‚
[ ]
ïÅ‚s2 śł
ðÅ‚ ûÅ‚
Przypadek 1 (¸ = 0.001, ¾ = 0.5, n = 500)
Przypadek 2 (¸ = 0.001, ¾ = 1.0 , n = 500)
Przypadek 3 (¸ = 0.001, ¾ = 6.3 , n = 2000)
Przypadek 4 (¸ = 0.001, ¾ = 6.0 , n = 2000)
Przypadek 5 (¸ = 0.01, ¾ = 6.0 , n = 200)
W przykładzie tym widzimy, jak przyjęcie zbyt du\
ego kroku czasowego ¸ powoduje na
¸
¸
¸
tyle znaczące pogorszenie dokładności rozwiązania numerycznego, \
e w praktyce nie
mogłoby być ono zaakceptowane.
Zadanie problemowe  drgania pręta
Jednorodny pręt AB o masie m jest obcią\
ony w środku geometrycznym dodatkową pionową
siłą Q = Q t . Końce A i B mogą przemieszczać się bez poślizgu odpowiednio wzdłu\
osi y i
( )
x. Punkty O, A oraz O, B połączone są odpowiednio pionową sprę\
yną o stałej ky i poziomą
sprÄ™\
yną o stałej kx, których długości w stanie naturalnym (nierozciągniętym) wynoszą
odpowiednio ly, lx. Przyjmując kąt Ć = Ć t jako współrzędną uogólnioną Lagrange a,
( )
wyprowadz
równanie ró\
niczkowe ruchu pręta z równania Lagrange a II-go rodzaju:
2
3ky
d Ć 3kx 3
= - îÅ‚
( ) ( ) ( )ûÅ‚ ( )-
( )
ðÅ‚l sin Ć - ly Å‚Å‚ cos Ć - ml îÅ‚lx - l cos Ć Å‚Å‚ sin Ć 2 ml P cos Ć (" )
ûÅ‚ ðÅ‚
dt2 ml
2 2
gdzie P = P t = G + Q t (G = m g , g oznacza przyśpieszenie ziemskie), l = lx + ly .
( ) ( )
Dopasowując się do oznaczeń argumentów i funkcji przyjętych w opisanych powy\
ej
2
dĆ d d
algorytmach: t = x , Ć = u , = v , ... = ... ', ... = ... " ostateczna postać
( ) ( ) ( ) ( )
dt dt dt2
równania ró\
niczkowego (" ) będzie zatem następująca:
3ky
3kx 3
u '' = - îÅ‚
( ) ( )-
( )ûÅ‚ ( ) ( )
ðÅ‚l sin u - ly Å‚Å‚ cos u ml îÅ‚lx - l cos u Å‚Å‚ sin u - 2 ml P cos u (" " )
ûÅ‚ ðÅ‚
ml
która z kolei jest równowa\
na układowi dwóch równań ró\
niczkowych zwyczajnych
u ' = v
v ' = h x,u,v
( )
gdzie
3ky
3kx 3
h x,u,v = - îÅ‚
( ) ( )- ( ) ( )ûÅ‚ ( ) ( )
ðÅ‚l sin u ly Å‚Å‚ cos u - ml îÅ‚lx - l cos u Å‚Å‚ sin u - 2 ml P cos u .
ûÅ‚ ðÅ‚
ml
Znajdz
numeryczne rozwiązanie równania ró\
niczkowego funkcjÄ… Runge_Kutta(...) dla
następujących danych:
x0 = 0 , u0 = 0.12435499454676195 , v0 = 0 , ¸ = 0.003 , n = 31999
f x,u,v = v , g x,u,v = h x,u,v
( ) ( ) ( )
m = 1.0 , kx = 2.0 , ky = 2.0 , lx = 4.0 , ly = 3.0 , g = 10.0
y
Ä™
ky
ly
P
kx
Ć ł
ź x
lx
Siłę Q = Q t = x zdefiniuj na cztery sposoby:
( )
I przypadek:
Q
t0
t
- G
II przypadek:
Q
G
t
t0
III przypadek: Q t = 0 = const
( )
IV przypadek: Q t = G sin 0.1Å"t
( ) ( )
gdzie t0 = 60 .
Wartości składowych wektorów x n +1 , U n +1 zapisz do pliku tekstowego i porównaj
[ ] [ ]
wczytane wartości z tymi podanymi na poni\
szych rysunkach.
Zale\
ność kąta Ć = Ć t kolejno dla I, II, III i IV przypadku
( )
ly
ëÅ‚ öÅ‚
Kąt Ć H" 0.643501 rad = arctan definiuje poło\
enie równowagi pręta dla przypadku
[ ]
ìÅ‚ ÷Å‚
lx
íÅ‚ Å‚Å‚
P = 0 (poło\
enie w którym sprę\
yny sÄ… w stanie naturalnym).
Kąt Ć H" -0.463648 rad definiuje poło\
enie równowagi pręta dla przypadku P = 2 G .
[ ]
Kąt Ć H" 0.124355 rad definiuje poło\
enie równowagi pręta dla przypadku P = G .
[ ]