Metoda Runge - Kutty

Równanie

 

t

,

x

f

x 



Rozwiązujemy stosując szereg Taylora





 

 

 

 























t

x

6

h

t

x

2

h

t

x

h

t

x

h

t

x

3

2

ale

   

t

,

x

f

t

x 



 

 

 





t

,

t

x

f

dt

d

x

dt

d

t

x









 





t

x

,

f

f

,

f

t

f

dt

dx

x

f

t

,

t

x

f

dt

d

















czyli

t

x

,

f

,

f

)

t

(

x







 

 

   

 





t

,

x

,

f

t

,

x

,

f

t

,

x

f

dt

d

x

dt

d

t

x

t

x













 





tt

xt

t

t

t

,

f

,

f

t

,

f

dt

dx

x

,

f

t

,

x

,

f

dt

d





















dt

,

df

f

,

f

dt

df

,

f

dt

d

x

x

x









dt

,

df

f

,

f

dt

df

,

f

dt

d

x

x

x





ale

t

x

,

f

f

,

f

dt

df





a

xt

xx

x

x

x

,

f

,

f

t

,

f

dt

dx

x

,

f

dt

,

df

















Podstawiając do

 

 

   

 





t

,

x

,

f

t

,

x

,

f

t

,

x

f

dt

d

x

dt

d

t

x

t

x













i porządkując mamy

tt

t

y

2

y

xt

xx

2

,

f

,

f

,

f

,

f

,

f

2

,

f

f

x















Metoda Runge -Kutty





p

2,3,...,

i

h

a

t

,

K

b

x

hf

K

t

,

x

hf

K

K

w

x

x

i

n

1

i

1

j

j

ij

n

i

n

n

i

p

1

i

i

i

n

1

n







































Sposób wyznaczania współczynników na przykładzie
metody drugiego rzędu (p=2):













h

b

t

,

t

,

x

hf

b

x

hf

w

t

,

x

hf

w

x

x

21

n

n

n

21

n

2

n

n

1

n

1

n















Drugi składnik rozwijamy w szereg Taylora
w otoczeniu punktu x

n

, t

n

















h

b

,

f

f

b

,

f

t

,

x

f

h

b

t

,

t

,

x

hf

b

x

f

21

t

21

x

n

n

21

n

n

n

21

n













Podstawiając i porządkując mamy:



 





 











n

n

t

n

n

n

n

x

2

21

2

n

n

2

1

n

1

n

t

,

x

,

f

t

,

x

f

t

,

x

,

f

h

b

w

t

,

x

hf

w

w

x

x















a porównując z szeregiem Taylora







 











n

n

t

n

n

x

n

n

2

n

n

n

1

n

t

,

x

,

f

t

,

x

,

f

t

,

x

f

2

h

t

,

x

hf

x

x











przy tych samych potęgach h otrzymujemy:

2

1

b

w

1

w

w

21

2

2

1







Przyjmując w

2

=1 mamy:

w

1

=0, b

21

=1/2

i stąd algorytm:

2

1

b

w

1

w

w

21

2

2

1















h

5

.

0

t

,

t

,

x

hf

5

.

0

x

hf

x

x

n

n

n

n

n

1

n











lub w1=w2=w i rozwiązując otrzymujemy:

w=0.5, b

21

=1 i stąd inny algorytm:

















1

n

n

n

n

n

n

n

1

n

t

,

t

,

x

hf

x

f

t

,

x

f

h

5

.

0

x

x













Przykład:
Dany jest dławik o charakterystce:

3

i

001

.

0

i

1

.

0 





Rezystancja dławika wynosi 0.5.

Obliczyć prąd płynący w obwodzie zasilanym
sem e(t)=100sin314t.

Schemat obwodu możemy przyjąć w postaci:

Suma spadków napięć pozwala zapisać równanie:

 

t

e

i

5

.

0

dt

d







Biorąc pod uwagę krzywą magnesowania:

dt

di

i

003

.

0

dt

di

1

.

0

dt

d

i

001

.

0

i

1

.

0

2

3













Podstawiając do równania obwodu i porządkując:

2

i

003

.

0

1

.

0

i

5

.

0

t

314

sin

100

dt

di







Warunek początkowy jest i

0

=i(t=0)=0.

Wybór kroku całkowania:

Stała czasowa liniowej części obwodu wynosi

0.1/0.5=0.2s.

Krok czasowy można przyjąć 0.2/10=20ms.

Okres wymuszenia T=20ms krok należy przyjąć
rzędu T/20=1ms. Prawdopodobnie będzie trzecia
harmoniczna więc przyjmujemy krok

h=0.2ms.

Obliczenia metodą Runge – Kutty według schematu:









2

n

1

n

1

n

2

n

n

1

K

x

x

h

5

.

0

t

,

K

5

.

0

x

hf

K

t

,

x

hf

K















 

2

i

003

.

0

1

.

0

i

5

.

0

t

314

sin

100

t

,

i

f







x=i;

Start: i(t=0)=i

0

=0

h=0.0002

i mamy:

0

i

003

.

0

1

.

0

i

5

.

0

0

314

sin

100

0002

.

0

K

2

0

0

1



























006279

.

0

K

K

5

.

0

i

003

.

0

1

.

0

K

5

.

0

i

5

.

0

0001

.

0

0

314

sin

100

0002

.

0

K

2

2

1

0

1

0

2















006279

.

0

K

i

i

2

0

1







t=h=0.0002

Metoda Runge – Kutty pozwala zmienić krok na
każdym etapie. Zwiększamy krok dwukrotnie.

h=0.0004

2

1

006279

.

0

003

.

0

1

.

0

006279

.

0

5

.

0

0002

.

0

314

sin

100

0004

.

0

K













02509

.

0

1

K 













2

1

1

2

K

5

.

0

006279

.

0

003

.

0

1

.

0

K

5

.

0

006279

.

0

5

.

0

0002

.

0

0002

.

0

314

sin

100

0004

.

0

K

















05003

.

0

K

2



i

2

=0.006279+K

2

i

2

=0.056312

Jak ocenić czy wolno zmienić długość kroku?
Czy zmniejszyć czy zwiększyć?

Ocena błędu metodą Rungego:

Dla metody rzędu p-go mamy:









 

2

p

1

p

n

n

1

n

1

n

h

O

h

t

,

x

x

t

x





















 

2

p

2

n

n

2

1

n

n

h

O

2

h

t

,

x

x

2

h

t

x



































 









 

2

p

1

p

n

n

2

1

2

n

1

n

h

O

2

h

t

,

x

2

x

t

x































stąd ocena błędu:





 

2

p

p

1

n

2

1

2

n

1

p

n

n

h

O

2

1

x

x

h

t

,

x























Znając ocenę błędu można poprawić rozwiązanie
podstawiając do









 

2

p

1

p

n

n

1

n

1

n

h

O

h

t

,

x

x

t

x

















 

2

p

p

1

n

2

1

2

n

p

1

n

p

1

n

h

O

1

2

x

x

2

x

x























lub dokładniej z równania:









 

2

p

1

p

n

n

2

1

2

n

1

n

h

O

2

h

t

,

x

2

x

t

x































 

2

p

p

1

n

2

1

2

n

2

1

2

n

dd

1

n

h

O

1

2

x

x

x

x

























W obliczanym przypadku musimy powtórzyć
obliczenia z krokiem 0.0002 i mamy dla t=0.0004:

K

1

=0.012545

K

2

=0.0188

i

1+1/2

=0.02508

Dla t=0.0006 mamy:

K

1

=0.025029

K

2

=0.031235

i

1+2*1/2

=0.056315

i

1+2*1/2

=0.056315

Obliczone z krokiem h=0.0004 było:

i

2

=0.056312

W tym przypadku p=2 i ze wzoru:





 

2

p

p

1

n

2

1

2

n

1

p

n

n

h

O

2

1

x

x

h

t

,

x























mamy oceną błędu:





00026

.

0

t

,

i

1

1





Rozwiązanie poprawione ze wzoru:

 

2

p

p

1

n

2

1

2

n

2

1

2

n

dd

1

n

h

O

1

2

x

x

x

x

























i

2

=0.0566316

Na wykonanie jednego kroku należało policzyć
funkcję f(i

n

,t

n

)

2 – h=0.0004
1+2 – h=0.0002
razem 5 - razy

Metoda IV –go rzędu





















4

3

2

1

n

1

n

n

3

n

4

n

2

n

3

n

1

n

2

n

n

1

K

K

2

K

2

K

6

1

x

x

h

t

,

K

x

hhf

K

h

5

.

0

t

,

K

5

.

0

x

hf

K

h

5

.

0

t

,

K

5

.

0

x

hf

K

t

,

x

hf

K

































Przy ocenie dokładności obliczeń metodą Rungego
wymaga 11-krotnego obliczenia f(x,t).

Metoda Mersona





 

30

K

K

8

K

9

K

2

h

error

h

t

,

2

K

4

K

3

K

x

hf

K

2

h

t

,

8

K

3

K

x

hf

K

3

h

t

,

6

K

K

x

hf

K

3

h

t

,

3

K

x

hf

K

t

,

x

hf

K

5

4

3

1

n

4

3

1

n

5

n

3

1

n

4

n

2

1

n

3

n

1

n

2

n

n

1





























































































6

K

K

4

K

x

x

5

4

1

n

1

n











tylko 5-cio krotne obliczanie f(x,t).

Przykład

Równanie wahadła:

0

sin

2











Niech =1s

-2

Warunki początkowe:





1

.

2

0

t









około 86°





0

0

t







Sprowadzamy do układu równań I-go rzędu















sin





Warunki początkowe:

0

1

.

2

0

0











Obliczenia chcemy prowadzić z dokładnością 0.001

Startujemy z krokiem h=0.1. Krok wybrano jako
0.1 okresu wahadła liniowego.

09972

.

0

3

K

sin

h

K

003324

.

0

3

K

h

K

09972

.

0

sin

h

K

0

h

K

1

0

2

1

0

2

0

1

0

1











































































099716

.

0

6

K

K

sin

h

K

003324

.

0

6

K

K

h

K

2

1

0

3

2

1

0

3































































099224

.

0

8

K

3

K

sin

h

K

0037394

.

0

8

K

3

K

h

K

3

1

0

4

3

1

0

4































































099701

.

2

K

4

K

3

K

sin

h

K

0098734

.

0

2

K

4

K

3

K

h

K

4

3

1

0

5

4

3

1

0

5







































































Błąd:

00013043

.

0

30

K

K

8

K

9

K

2

00032914

.

0

30

K

K

8

K

9

K

2

5

4

3

1

1

5

4

3

1

1





















































Dokładność założona została osiągnięta.
W następnym kroku można zwiększyć krok.

Rozwiązanie w chwili t=0.1

099386

.

6

K

K

4

K

49186

.

1

6

K

K

4

K

5

4

1

0

1

5

4

1

0

1











































i do następnego kroku możemy wystartować z nową
wartością kroku h

Metody włożone

lub

Metody Fehelberga – Runge -Kutty

Stosujemy metodę Runge – Kutty rzędu p i rzędu p+1
i aby zmniejszyć liczbę obliczanych współczynników
wybieramy je tak, że w obu metodach jest pierwszych
p współczynników K jednakowe, czyli





n

n

1

t

,

x

hf

K 

























h

c

t

,

K

a

x

hf

K

i

n

1

i

1

j

j

ij

n

i

i=2,3,..,p+1

i mamy dla metody rzędu p-go











p

1

i

i

p

i

n

1

n

K

w

x

x

a dla metody rzędu (p+1)-go















1

p

1

i

i

1

p

i

n

1

n

K

w

x

x

Ocenę błędu można zrobić stosunkowo prosto





 

2

p

1

p

n

n

p

p

1

i

i

p

i

n

1

n

h

O

h

t

,

x

K

w

x

x

























 

3

p

2

p

n

n

1

p

1

p

1

i

i

1

p

i

n

1

n

h

O

h

t

,

x

K

w

x

x



























Po odjęciu stronami otrzymujemy:









 

2

p

1

p

1

i

i

p

i

1

p

i

n

n

p

h

O

K

w

w

t

,

x



















gdzie

0

w

p

1

p





Znając błąd możemy postępować jak w metodzie
Mersona i rozwiązanie przyjmować z dokładniejszej
metody rzędu p+1.

Najczęściej stosowana metoda RKF45 ma współczynniki



























































h

13

12

t

,

2197

K

7296

K

7200

K

1932

x

hf

K

h

8

3

t

,

32

K

9

K

3

x

hf

K

h

25

.

0

t

,

K

25

.

0

x

hf

K

t

,

x

hf

K

n

3

2

1

n

4

n

2

1

n

3

n

1

n

2

n

n

1



















































2

h

t

,

K

40

11

K

4104

1859

K

2565

3544

K

2

K

27

8

x

hf

K

h

t

,

K

4104

845

K

513

3680

K

8

K

216

439

x

hf

K

n

5

4

3

2

1

n

6

n

4

3

2

1

n

5

6

5

4

3

1

n

K

55

2

K

50

9

5

1

K

4104

2197

56430

28561

K

2565

1408

12825

6656

K

216

25

135

16













 























































Błąd

Rozwiązanie wykorzystując metodę dokładniejszą jest

6

5

3

2

1

n

1

n

K

55

2

K

50

9

K

56430

28561

K

12825

6656

K

135

16

x

x

















Metoda gwarantuje obliczenia z błędem rzędu h

4

.

Przykład

Metody Rungego – Kutty a równania sztywne

0

x

x

01

.

1

x











z warunkiem początkowym:

 

 

0

0

x

1

0

x







Rozwiązanie analityczne jest:

 





 

t

exp

99

1

t

01

.

0

exp

99

100

t

x









i wykres jest:

0

1.5

3

4.5

6

0.96

0.968

0.976

0.984

0.992

1

x t

( )

t

dla t[0,5]

0

150

300

450

600

0

0.2

0.4

0.6

0.8

1

x 

 



dla t[0,600]

Wyniki otrzymane z metody Rungego – Kutty IV rzędu

dla t[0,10] i h=0.05

0

2

4

6

8

10

0.9

0.93

0.95

0.98

1

w

n

x n h



(

)

n h



i błąd

 

 

 

100

nh

x

nh

x

w

n

bl

n





0 40 80 120 160 200

0

610

9

1.210

8

1.810

8

2.410

8

310

8

bl n

( )

n

Po czasie t=10s można pominąć

 





 

t

exp

99

1

t

01

.

0

exp

99

100

t

x









drugi wyraz w rozwiązaniu czyli praktycznie

 





t

01

.

0

exp

99

100

t

x





a więc wystarczy liczyć do 1000s z krokiem równym 100/20=5s

i mamy:

0

2

4

6

8

10

1.510

5

1.110

5

710

4

310

4

110

4

510

4

v

n

n

czyli rozwiązanie się rozbiega. Po 10 krokach czyli dla

t=10+5·10=60s wynik jest z błędem 10

7

%

dla kroku dwa razy mniejszego czyli h=2.5s mamy:

0

200 400 600 800 1000

0

0.2

0.4

0.6

0.8

1

v

n

x 10 n hn





(

)

n hn



bardzo dobre rozwiązanie. Błąd jest:

 









100

nh

10

x

nh

10

x

v

n

bl

n









0

80 160 240 320 400

0

110

5

210

5

310

5

410

5

b n

( )

n

Błąd procentowy jest do przyjęcia. Otrzymany wynik pokazuje,

że istnieje granica stabilności absolutnej dla metod typu

Rungego - Kutty

Wykazano, że warunek stabilności absolutnej dla metody IV
rzędu jest:

78

.

2

h

min





gdzie 

min

jest najmniejszą stałą czasową w układzie równań

sztywnych

W analizowanym przypadku:

 





 

t

exp

99

1

t

01

.

0

exp

99

100

t

x









stałe czasowe są: 

min

=1s i 

max

=100s, a więc dopuszczalny krok

jest h

dop

=2.78. Dla kroku h=5>h

dop

mamy utratę stabilności,

a dla h=2.5<h

dop

mamy prawidłowy przebieg obliczeń.