Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-1

Równania różniczkowe zwyczajne
Zagadnienie początkowe

Szukamy różniczkowalnej y(x) spełniającej powyższe warunki.
Warunek Lipschitza:

2

1

2

1

2

1

)

,

(

)

,

(

0

,

],

,

[

y

y

L

x

y

f

x

y

f

L

R

y

y

b

a

x

−

≤

−

>

∃

∈

∈

∀

Jeżeli ciągła funkcja

R

R

b

a

f

→

×

]

,

[

:

spełnia warunek Lipschitza:

to zagadnienie początkowe ma dokładnie jedno różniczkowalne w
sposób ciągły rozwiązanie y:[a,b]→R.

c

a

y

x

x

y

f

dx

dy

=

=

)

(

),

),

(

(

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-2

SCHEMAT
RÓŻNICOWY

rozwiązanie numeryczne: na przedziale [a,b]
punkty:

a=x

0

, x

1

, x

2

, ....., x

n

=b,

x

i

-x

i-1

=h

i

wartości przybliżone: y

0

, y

1

, y

2

, ....., y

n

wartości dokładne: y(x

0

), y(x

1

), .......,y(x

n

)

Błędy

SCHEMAT
RÓŻNICOWY

ANALIZA BŁĘDU

?

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-3

Wybrane schematy różnicowe jednokrokowe

Metoda EULERA

)

x

,

y

(

f

h

y

y

)

x

,

y

(

f

h

y

y

),

x

),

x

(

y

(

f

dx

dy

n

n

n

n

n

n

n

n

1

1

1

1

+

+

+

+

=

−

−

≈

−

=

h

)

x

,

y

(

f

y

y

h

)

x

,

y

(

f

y

y

n

n

n

n

n

n

n

n

1

1

1

1

+

+

+

+

+

=

+

=

otwarta

(jawna) zamknięta (niejawna)

Zmodyfikowana metoda Eulera
pół kroku w kierunku pochodnej - poprawienie pochodnej - cały krok w
poprawionym kierunku

)

h

x

,

y

(

f

f

h

f

y

y

)

x

,

y

(

f

f

n

n

n

n

n

n

n

n

n

2

2

2

1

2

1

2

1

+

=

+

=

=

+

+

+

h

f

y

y

n

n

n

2

1

1

+

+

+

=

Schemat jednokrokowy

h

)

h

,

x

,

y

,

y

(

y

y

n

n

n

f

n

n

1

1

+

+

Φ

+

=

- niejawny (zamknięty)

h

)

h

,

x

,

y

(

y

y

n

n

f

n

n

Φ

+

=

+1

- jawny (otwarty)

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-4

Schematy Rungego-Kutty r-poziomowy (etapowy) schemat R-K

)

h

,

x

,

y

(

K

c

)

h

,

x

,

y

(

n

n

i

r

i

i

n

n

f

∑

=

=

Φ

1

r

,...,

,

i

),

b

h

x

,

K

b

h

y

(

f

)

h

,

x

,

y

(

K

r

j

ij

j

r

j

ij

i

2

1

1

1

=

+

+

=

∑

∑

=

=

otwarty (jawny):

i

j

b

ij

≥

= 0

r

,...,

i

),

b

h

x

,

K

b

h

y

(

f

)

h

,

x

,

y

(

K

),

x

,

y

(

f

K

i

j

ij

j

i

j

ij

i

2

1

1

1

1

1

=

+

+

=

=

∑

∑

−

=

−

=

RK4:

)

K

K

K

K

(

h

y

y

),

h

x

,

K

y

(

f

K

),

h

x

,

K

y

(

f

K

),

h

x

,

K

y

(

f

K

),

x

,

y

(

f

K

n

n

n

n

n

n

n

n

n

n

4

3

2

1

1

3

4

2

3

1

2

1

2

2

6

1

2

1

2

1

2

1

2

1

+

+

+

+

=

+

+

=

+

+

=

+

+

=

=

+

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-5

Zbieżność metod jednokrokowych i sterowanie długością kroku

h

)

h

,

x

,

y

(

y

y

n

n

f

n

n

Φ

+

=

+1

h

)

h

,

x

),

x

(

y

(

)

x

(

y

y

n

n

f

n

n

Φ

+

=

+1

]

h

)

h

,

x

),

x

(

y

(

)

x

(

y

[

)

h

x

(

y

r

n

n

f

n

n

n

Φ

+

−

+

=

+1

- błąd lokalny

Metoda rządu p:

)

h

(

O

h

)

x

),

x

(

y

(

)

h

(

r

p

p

n

n

n

2

1

1

+

+

+

+

=

ϕ

(

1

1

)

(

+

+

<

p

n

Ch

h

r

)

metoda p
Eulera 1
zmodyfikowana Eulera

2

RK2,3,4 2,3,4
RKm m=5,6,..

p<m

Jeżeli metoda jednokrokowa jest rzędu p i spełnia warunek Lipschitza

y

y

L

h

x

y

h

x

y

f

f

~

)

,

,

~

(

)

,

,

(

−

≤

Φ

−

Φ

Φ

to błąd globalny można oszacować przez

(

)

l

l

n

l

p

l

l

n

l

x

x

h

Ch

x

y

y

−

=

≤

−

+

−

=

=

1

1

,

,

1

,

0

max

max

,

,

1

,

0

max

,

)

(

max

"

"

.

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-6

1

+

→

n

h

n

n

y

)

x

(

y

,

x

2

1

2

1

2

2

+

+

→

→

n

h

h

n

n

y

)

x

(

y

,

x

...

h

y

)

h

x

(

y

p

n

n

+

=

−

+

+

+

1

1

ϕ

...

h

y

)

h

x

(

y

p

n

n

+

⎟

⎠

⎞

⎜

⎝

⎛

=

−

+

+

+

+

1

2

1

2

1

2

2

ϕ

ERR=

)

y

y

(

h

n

n

p

p

p

1

2

1

2

1

1

1

2

2

+

+

+

+

−

−

≈

ϕ

dla p=4 : ERR=

)

y

y

(

n

n

1

2

1

2

1

15

16

+

+

+

−

ERR < y

max

RELREER+ABSERR

15

16

1

2

1

2

1

1

+

+

+

+

−

=

⇒

n

n

n

y

y

:

y

ERR > y

max

RELREER+ABSERR

⇒

zmniejszyć h

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-7

inaczej

ERR <

2

1

+

+

n

n

y

y

RELREER+ABSERR lub

ETOL=

1

2

1

<

+

+

+

ABSERR

RELERR

y

y

ERR

n

n

szukamy

α

, by

1

≈

)

h

(

ETOL

α

:

)

h

(

ETOL

)

h

(

ETOL

5

α

α

≈

5

1

)

h

(

ETOL

=

α

0.9

,

h

h

nowe

α

=

Jeżeli 1

<

)

h

(

ETOL

, to

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

5

9

0

5

,

)

h

(

ETOL

.

min

α

Metoda wymaga dla RK4 4+7=11 obliczeń prawej strony.

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-8

Algorytm Rungego-Kutty-Fehlberga stosuje dwa schematy
Rungego-Kutty m+1 i m etapowy z odpowiednio dobranymi
współczynnikami. Schemat m-etapowy jest rzędu p a schemat m+1
etapowy jest rzędu p+1, a współczynniki

i

K są jednakowe.

h

)

h

,

x

,

y

(

y

y

n

n

f

n

n

Φ

+

=

+1

=

)

h

,

x

,

y

(

K

c

y

n

n

i

m

i

i

n

∑

+

=

+

1

1

h

)

h

,

x

,

y

(

y

y

~

n

n

f

n

n

Φ

+

=

+1

=

)

h

,

x

,

y

(

K

c~

y

n

n

i

m

i

i

n

∑

=

+

1

0

1

=

+

m

c~

Wtedy

)

c~

c

(

h

)

h

(

ERR

i

m

i

i

−

=

∑

+

=

1

1

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-9

Metody wielokrokowe

Metody Adamsa: Adamsa-Bashfortha (jawne), Adamsa –Moultona
(niejawne)

∫

+

+

=

+

1

1

n

n

x

x

n

n

dx

)

x

),

x

(

y

(

f

)

x

(

y

)

x

(

y

przybliżamy funkcję podcałkową wielomianem interpolacyjnym
Lagrange’a w węzłach

⎩

⎨

⎧

=

−

=

−

−

−

jawna

metoda

k

,...

,

i

niejawna

metoda

k

,...

,

,

i

))

x

,

y

(

f

,

x

(

i

n

i

n

i

n

1

0

1

0

1

całkujemy wielomian i otrzymujemy wzory postaci

)

x

,

y

(

f

~

h

y

~

y

k

j

j

n

j

n

j

j

n

k

j

j

n

∑

∑

=

−

−

−

=

+

+

=

0

0

1

β

α

dla metody jawnej

)

x

,

y

(

f

h

)

x

,

y

(

f

h

y

y

n

n

k

j

j

n

j

n

j

j

n

k

j

j

n

1

1

1

0

0

1

+

+

−

=

−

−

−

=

+

+

+

=

∑

∑

β

β

α

dla

metody niejawnej

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-10

Metoda wstecznego różniczkowania
przybliżamy rozwiązanie y(x) wielomianem interpolacyjnym W(x)
zbudowanym na węzłach

⎩

⎨

⎧

=

−

=

−

−

jawna

metoda

k

,...

,

i

niejawna

metoda

k

,...

,

,

i

)

y

,

x

(

i

n

i

n

1

0

1

0

1

obliczamy pochodną tego wielomianu W’(x) , która przybliża
pochodną rozwiązania,
z równości

)

(

'

)

,

(

)

(

n

n

n

n

x

W

x

y

f

x

y

≈

=

′

otrzymujemy wzór

)

,

(

~

~

0

0

1

n

n

j

n

k

j

j

n

x

y

f

h

y

y

β

α

+

=

−

=

+

∑

dla metody jawnej

a z

)

(

'

)

,

(

)

(

1

1

1

1

+

+

+

+

≈

=

′

n

n

n

n

x

W

x

y

f

x

y

)

,

(

1

1

1

0

1

+

+

−

−

=

+

+

=

∑

n

n

j

n

k

j

j

n

x

y

f

h

y

y

β

α

dla metody niejawnej

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-11

Stosowanie metod niejawnych

)

x

,

y

(

f

h

)

x

,

y

(

f

h

y

y

n

n

k

j

j

n

j

n

j

j

n

k

j

j

n

1

1

1

0

0

1

+

+

−

=

−

−

−

=

+

+

+

=

∑

∑

β

β

α

jest nieliniowym równaniem algebraicznym względem

1

+

n

y

i jest

rozwiązywane metodą iteracyjną (najczęściej iteracji prostej):

)

x

,

y

(

f

h

)

x

,

y

(

f

h

y

y

n

]

i

[

n

k

j

j

n

j

n

j

j

n

k

j

j

]

i

[

n

1

1

1

1

0

0

1

+

−

+

−

=

−

−

−

=

+

+

+

=

∑

∑

β

β

α

, i=1,2,...

Wymaga ona punktu startowego, np. z metody Eulera

)

x

,

y

(

hf

y

y

n

n

n

]

[

n

+

=

+

0

1

. Uzyskanie dokładnego rozwiązania wymaga

wtedy wielu iteracji.
Metody typu PREDYKTOR-KOREKTOR wyznaczają przybliżenie
początkowe jawną metodą wielokrokową o takiej samej liczbie
kroków, a następnie stosują kilka iteracji rozwiązujących równanie
nieliniowe.
Algorytm GEAR’A

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-12

Liniowe układy równań różniczkowych

0

0

y

)

(

y

),

x

(

y

A

)

x

(

y

dx

d

G

G

G

G

=

=

gdzie A jest macierzą nxn o różnych wartościach własnych

n

,...,

,

i

i

2

1

=

λ

i wektorach własnych

n

,...,

,

i

v

i

2

1

=

G

Wartość własna i wektor własny:

0

≠

=

i

i

i

i

v

,

v

v

A

G

G

G

λ

(

)

0

=

−

i

i

v

A

I

G

λ

(

)

0

=

− A

I

det

i

λ

czyli wartości własne są pierwiastkami równania

(

)

0

=

− A

Is

det

wielomian stopnia n (wielomian charakterystyczny A)

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-13

Niech:

[

]

n

v

v

v

V

G

"

G

G

2

1

=

,

)

x

(

y

)

x

(

z

V

G

G

=

Wtedy:

Λ

= V

AV

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

=

Λ

n

λ

λ

λ

0

0

0

0

0

0

2

1

"

#

%

#

#

"

"

Λ

=

−

AV

V

1

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-14

Niech:

)

x

(

y

)

x

(

z

V

G

G

=

0

0

1

1

0

0

z

y

V

)

(

y

V

)

(

z

),

x

(

z

AV

)

x

(

z

dx

d

V

G

G

G

G

G

G

=

=

=

=

−

−

0

0

1

1

1

0

0

z

y

V

)

(

y

V

)

(

z

),

x

(

z

AV

V

)

x

(

z

dx

d

G

G

G

G

G

G

=

=

=

=

−

−

−

0

0

1

1

2

1

0

0

0

0

0

0

0

0

z

y

V

)

(

y

V

)

(

z

),

x

(

z

)

x

(

z

dx

d

n

G

G

G

G

G

"

#

%

#

#

"

"

G

=

=

=

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

=

−

−

λ

λ

λ

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-15

n

,...,

,

i

z

)

(

z

),

x

(

z

)

x

(

z

dx

d

i

i

i

i

i

2

1

0

0

=

=

=

λ

n

,...,

,

i

z

e

)

x

(

z

i

x

i

i

2

1

0

=

=

λ

),

(

z

e

e

e

)

x

(

z

x

x

x

n

0

0

0

0

0

0

0

2

1

G

"

#

%

#

#

"

"

G

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

λ

λ

λ

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-16

),

(

z

e

e

e

V

)

x

(

z

V

x

x

x

n

0

0

0

0

0

0

0

2

1

G

"

#

%

#

#

"

"

G

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

λ

λ

λ

),

(

y

V

e

e

e

V

)

x

(

y

x

x

x

n

0

0

0

0

0

0

0

1

2

1

G

"

#

%

#

#

"

"

G

−

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

λ

λ

λ

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

−

T

n

T

T

w

w

w

:

V

G

#

G

G

2

1

1

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-17

),

(

y

w

v

e

)

x

(

y

n

i

T

i

i

x

i

0

1

G

G

G

G

∑

=

=

λ

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-18

Stabilność i sztywność

Rozważmy układ n równań

0

0

y

)

(

y

),

x

(

Ay

dx

dy

=

=

rozwiązanie:

∑

=

=

=

n

i

x

i

Ax

i

e

c

y

e

)

x

(

y

1

0

λ

dąży do 0 dla

0

<

i

Re

λ

Metodą Eulera

n

n

n

n

y

)

hA

I

(

hAy

y

y

+

=

+

=

+1

:

0

1

y

)

hA

I

(

y

+

=

,

0

2

1

2

y

)

hA

I

(

y

)

hA

I

(

y

+

=

+

=

0

3

2

3

y

)

hA

I

(

y

)

hA

I

(

y

+

=

+

=

.............................................

n

,...,

i

h

y

lim

i

n

n

1

1

1

0

=

<

+

⇔

=

∞

→

λ

Obszar stabilności schematu różnicowego: zbiór wartości
zespolonych

h

i

λ

, dla których wszystkie rozwiązania zadania

testowego są ograniczone dla

∞

→

n

. Jeżeli obszar stabilności

zawiera punkt 0, to metodę nazywamy stabilną.
Metoda A-stabilna A(a)-stabilna , A(

α) stabilna

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-19

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-20

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-21

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-22

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-23

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-24

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-25

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-26

Obszary stabilności absolutnej

Adamsa-Bashfortha
l. kroków

1

2

3

4

rząd 1

2

3

4

q

*

-2 -1 -6/11

-3/10

Adamsa-Moultona
l. kroków

1

2

3

4

rząd 2

3

4

5

q

*

-

∞

-6 -3 -90/49

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-27

Wstecznego różniczkowania
k 1

2

3

4

5

6

a 0

0

-0.1

-0.7

-2.4 -6.1

α

90 90 88 73 51 18

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-28

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-29

Rungego-Kutty
p=m 1

2

3 4

Lewa
granica
obszaru
stabilności

-2 -2 -2.51 -2.78

Żadna metoda jawna nie jest A-
stabilna A(a)-stabilna lub A(α)-
stabilna dla żadnych a i α

.

Każda metoda jawna ma ograniczony obszar stabilności.

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-30

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-31

ode45 is based on an explicit Runge-Kutta (4,5) formula, the
Dormand-Prince pair. In general, ode45 is the best solver to apply
as a first try for most problems. For this reason, ode45 is the default
solver used by Simulink for models with continuous states.

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-32

ode23 is also based on an explicit Runge-Kutta (2,3) pair of Bogacki
and Shampine. It can be more efficient than ode45 at crude
tolerances and in the presence of mild stiffness.
ode113 is a variable-order Adams-Bashforth-Moulton PECE solver.
It can be more efficient than ode45 at stringent tolerances.
ode15s is a variable order solver based on the numerical
differentiation formulas (NDFs). These are related to but are more
efficient than the backward differentiation formulas, BDFs. If you
suspect that a problem is stiff, or if ode45 failed or was very
inefficient, try ode15s.
ode23s is based on a modified Rosenbrock formula of order 2.
Because it is a one-step solver, it can be more efficient than ode15s
at crude tolerances. It can solve some kinds of stiff problems for
which ode15s is not effective.
ode23t is an implementation of the trapezoidal rule using a "free"
interpolant. Use this solver if the problem is only moderately stiff
and you need a solution without numerical damping.

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-33

ode23tb is an implementation of TR-BDF2, an implicit Runge-Kutta
formula with a first stage that is a trapezoidal rule step and a
second stage that is a backward differentiation formula of order
two. By construction, the same iteration matrix is used in evaluating
both stages. Like ode23s, this solver can be more efficient than
ode15s at crude tolerances.

Note: For a stiff problem, solutions can change on a time scale that is
very short compared to the interval of integration, but the solution of
interest changes on a much longer time scale. Methods not designed for
stiff problems are ineffective on intervals where the solution changes
slowly because they use time steps small enough to resolve the fastest
possible change. Jacobian matrices are generated numerically for ode15s
and ode23s. For more information, see Shampine, L. F., Numerical
Solution of Ordinary Differential Equations, Chapman & Hall, 1994.

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-34

Przykłady
Rozwiązanie wolnozmienne

10

2

1

3

)

(

,

50

0

,

1

)

0

(

3

)

(

1

)

(

1

.

0

)

(

t

e

t

f

t

f

t

f

t

f

dt

t

df

−

+

=

<

<

=

⎟

⎠

⎞

⎜

⎝

⎛ −

=

10

-12

10

-10

10

-8

10

-6

10

-4

10

-2

10

-14

10

-12

10

-10

10

-8

10

-6

10

-4

10

-2

Błąd w funkcji RelTol

RKF23
RKF45
Adam s PC
Gear W R

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-35

10

-12

10

-10

10

-8

10

-6

10

-4

10

-2

10

1

10

2

10

3

10

4

Kroki w funkcji RelTol

RKF23
RKF45
Adam s PC
Gear W R

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 9

W9-36

10

-12

10

-10

10

-8

10

-6

10

-4

10

-2

10

1

10

2

10

3

10

4

10

5

LW F w funkcji RelTol

RKF23
RKF45
Adam s PC
Gear W R

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne ET3 wykład 6

W6-37

Równanie sztywne

(

)

)

(

)

(

,

2

1

3

)

(

,

50

0

500

,

2

)

0

(

)

(

)

(

)

(

)

(

10

t

f

e

t

y

e

t

f

t

a

y

dt

t

df

t

f

t

y

a

dt

t

dy

at

t

+

=

+

=

<

<

=

=

+

−

−

=

−

−

0

5

10

15

20

25

30

35

40

45

50

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne ET3 wykład 6

W6-38

10

-6

10

-5

10

-4

10

-3

10

-2

10

1

10

2

10

3

10

4

10

5

LW F w funkcji RelTol

RKF23
RKF45
Adams PC
Gear W R