Szereg zespolony.

Dana jest funkcja określona przez podanie jej wartości f

n

w punktach:

N

n

2

x

n





gdzie n=0,1,2,...,N-1. Aproksymujemy funkcję wielomianem
trygonometrycznym postaci:

1

i

gdzie

)

ikx

exp(

c

)

x

(

w

1

N

0

k

k













Otrzymujemy N równań dla wyznaczenia współczynników c

k

:











1

N

0

k

n

k

n

n

)

ikx

exp(

c

f

)

x

(

w

Rozwiązanie ostatniego układu równań czyli współczynniki c

k

są określane równaniami:











1

N

0

n

n

n

p

)

ipx

exp(

f

N

1

c

Idea szybkie transformaty Fouriera tzw. FFT

Fast Fourier Transform

Ponieważ

N

n

2

x

n





więc

























N

pn

2

i

exp

ipx

exp

n

Oznaczmy:



















N

i

2

epx

w

Zauważmy, że

Re

Im

w=w

N+1

N

2

w

p

=w

p+N

Każda całkowita potęga w leży na okręgu jednostkowym
i co więcej jeżeli wykładnik p potęgi w

p

jest większy od N

to punkty się nakrywają. Na tym spostrzeżeniu bazuje FFT.

Piszemy:









1

N

0

n

n

pn

p

f

w

N

1

c

Możemy zapisać w postaci macierzowej:

Oznaczając:

1

N

,...,

1

,

0

n

dla

N

f

g

n

n







 

 

 

n

pn

p

g

w

c 

gdzie

 





































2

1

N

1

N

0

1

N

1

0

0

0

0

pn

w

w

w

.

.

.

w

w

w

w

w

w

w









Ponieważ w

0

=1 więc nie będziemy pisać zerowej kolumny i wiersza.

Dalej mamy związki:



















N

i

2

epx

w





 

k

k

N

w

k

N

i

2

exp

k

N

i

2

exp

N

N

i

2

exp

k

N

N

i

2

exp

w

















 













 



































czyli

 

k

k

N

w

w







a więc wiersze i kolumny: 1 i N-1
2 i N-2
..............
k i N-k
..............
N/2-1 i N/2+1 dla N parzystych
(N-1)2 i (N+1)/2 dla N nieparzystych

są sprzężone

.

W praktyce najczęściej stosowane N=2

M

.

Jeżeli liczba węzłów interpolacyjnych mniejsza od 2

M

,

to uzupełniamy zerami.

N=8

w

w

2

w

3

w

4

=-1 (w

*

)

3

(w

*

)

2

w

*

w

2

w

4

w

6

w

8

=1 (w

*

)

6

(w

*

)

4

(w

*

)

2

w

3

w

6

w

9

w

12

=-

1

(w

*

)

9

(w

*

)

6

(w

*

)

3

w

4

w

8

w

12

w

16

=1 (w

*

)

1

2

(w

*

)

8

(w

*

)

4

w

5

w

10

w

15

w

20

=-

1

(w

*

)

1

5

(w

*

)

1

0

(w

*

)

5

w

6

w

12

w

18

w

24

=1 (w

*

)

1

8

(w

*

)

1

2

(w

*

)

6

w

7

w

14

w

21

w

28

=-

1

(w

*

)

2

1

(w

*

)

1

4

(w

*

)

7



















8

i

2

epx

w











 



8

i

2

epx

w

*

lub inaczej

a b c

-

1

c

*

b

*

a

*

b -

1

b

*

1 b -1 b

*

c b

*

a

-

1

a

*

b c

*

-1 1 -1 1 -1 1 -1

c

*

b a

*

-

1

a b

*

c

b

*

-

1

b 1 b

*

-1 b

a

*

b

*

c

*

-

1

c

b a

















































7

6

5

4

3

2

1

0

f

f

f

f

f

f

f

f

dla otrzymania tablicy mnożników wystarczy obliczyć połowę

pierwszego wiersza!!!

np. a=a

r

+ia

i

oraz a

*

=a

r

-ia

i

czyli np.

af

1

+a

*

f

7

=a

r

(f

1

+f

7

)+ia

i

(f

1

-f

7

)

i podobnie dla innych

operacji.

Wykorzystanie przedstawionych uproszczeń pozwala w stosunku

do zwykłego algorytmu zawierającego N

2

działań zespolonych

zmniejszyć ich liczbę dla N=2

M

do 2NM

1

2

3

4

5

6

0

1200

2400

3600

4800

6000

2

2 m



m2

m 1





m

Rozwiązywanie równań algebraicznych

f(x)=0

Metoda bisekcji

Przykład:

0

1

x

4

x

3







x

-1 0

-0.5

-0.25

-0.125

-

0.1875

f(x) -4 1

-

1.125

-

0.01562

5

0.4980

45

0.2430

8

x

-

0.21875

-

0.23437

5

-

0.242187

5

-

0.24609375

f(x

)

0.11453

2

0.04962
54

0.017044
5

0.00072103

7

x

-
0.248046
9

-
0.2465820
3

-
0.2463379

-
0.2462158

f(x) -

0.007449
1

-
0.0013209
8

-
0.0002999
3

0.0002105
7

Zaleta metody: Jeżeli pierwiastek istnieje, to go znajdziemy.

Wada metody: Duża liczba obliczeń

Regula falsi.

Założenia:

a) funkcja ma w przedziale [a,b] tylko jeden pierwiastek
i zachodzi f(a)f(b)<0,
b) jest funkcja jest klasy C

2

[a,b], pierwsza i druga pochodna

nie zmieniają znaku na przedziale [a,b].

Funkcja spełniająca powyższe założenia musi mieć w otoczeniu

miejsca zerowego jeden z następujących przebiegów:

f(a)

a

b

f(b)

x

y

f(a)

a

b

f(b)

x

y

f(a)

a

b

f(b)

x

y

f(a)

a

b

f(b)

x

y

Przebieg obliczeń metodą regula falsi:

x

y

a

b

f(a)

f(b)

x

1

f(x

1

)

x

2

analitycznie: ustalamy koniec z
warunku f(x

1

)f(a)<0 lub f(x

1

)f(b)<0

Prowadzimy prostą:

   

   

a

f

b

f

a

f

x

f

a

b

a

x











   

   

a

f

b

f

a

f

x

f

a

b

a

x











ale f(x

1

)=0 stąd

 

   

a

f

b

f

a

f

a

b

a

x

1











lub

     

a

f

a

f

b

f

a

b

a

x

1









Dla n-tej iteracji mamy b=x

n-1

i podstawiając mamy:



    

a

f

a

f

x

f

a

x

a

x

1

n

1

n

n













Ocena błędu dla dostatecznie małego przedziału [x

n-1

,x

n

]

można przyjąć jako:

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

d













Metoda regula falsi jest zbieżna dowolnej funkcji ciągłej
na przedziale [a,b].
Poszukiwanie pierwiastka zostaje zakończone jeżeli:







 1

n

n

x

x

Metoda jest wolno zbieżna.

Przykład:

0

1

x

4

x

3







x

-1 0

-0.2

-

0.2366412

-

0.24423354

f(x) -4 1 0.19

2

0.0401342

7

0.00849729

2

     

 

   

2

.

0

x

4

4

1

1

0

1

x

a

f

a

f

b

f

a

b

a

x

1

1

1





























Ponieważ f(-1)=-4, a f(x

1

)=0.192,

więc stałym punktem będzie x=-1

x

-

0.24583563

2

-

0.2461749

2

-

0.24624682

9

f(x) 0.00180035

9

0.0003816

03

0.00008089

1

x

-

0.24626207

2

f(x) 0.00001714

7

w metodzie bisekcji potrzebowaliśmy

14 kroków

ocena błędu:

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

d













)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

d













ocena błędu:





6

d

10

1

.

4

246262072

.

0

x











Metoda siecznych

Przepis:

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

1

n















Przykład:

0

1

x

4

x

3







x

-1

0

-0.2

-

0.24752475

2

-

0.24625643

9

f(x) -4

1 0.192

-

0.00526448

1

0.00004070

3

w regula falsi potrzeba 8 kroków

x

-

0.2462661

7

f(x)

0.907E-8

w 6-tym kroku

Koniecznie trzeba obliczać f(x

n

) i jeżeli zaczyna narastać

należy zawęzić przedział i powtórzyć obliczenia.

Niebezpieczeństwo znalezienia fałszywego pierwiastka.

Metoda szybsza niż reguła falsi.

a

b

x

1

Pierwsza iteracja musi startować

z punktów spełniających warunek:

f(a)f(b)<0

Metoda Newtona - Raphsona

Niech  małe w mamy:



  

 

 

2

x

f

x

f

x

f

x

f

2



















Pomijając małe drugiego rzędu 

2

mamy, że f(x+)=0,

jeżeli

 

 

x

f

x

f









Graficznie:

x

y

x

n



n

 

n

n

x

f

tg







Równanie prostej stycznej
w punkcie x

n

jest:

  

  

n

n

n

x

f

x

x

x

f

y









x

n+1

Prosta:

  

  

n

n

n

x

f

x

x

x

f

y









przechodzi przez zero, czyli y=0, w punkcie x

n+1

i mamy:

 

 

n

n

n

1

n

x

f

x

f

x

x









Przykład:

0

1

x

4

x

3







 

x

0

-0.25

-
0.24626865
7

-
0.24626617
2

f(x) 1 -

0.01562
5

-
0.00001039

-0.2E-10

W 3 krokach dokładność osiągana w metodzie siecznych

w 5 krokach

W obliczeniach numerycznych pochodną najczęściej oblicza się
numerycznie:

Metoda Newtona – Raphsona jest zbieżna kwadratowo, tzn.

 

 

i

i

2
i

1

i

x

f

2

x

f















  

  

h

x

f

h

x

f

x

f

i

i

i









„Pechowe” przypadki:

x

f(x)

x

0

x

1

x

2

rozbieżna

Zmniejszyć przedział

[x

d

,x

0

]

x

d

cykl

x

f(x)

x

1

=x

3

=...

x

2

=x

4

=...

x

d

Budując procedurę należy się zabezpieczyć przed taką możliwością.

Wystartować z punktu x

1

znajdującego się bliżej x

d

Pierwiastki wielokrotne:

 

 

0

x

f

i

0

x

f

d

d







Przy pierwiastkach wielokrotnych badać funkcję:

)

x

(

f

)

x

(

f

)

x

(

g





Pierwiastki zespolone

Przykład

0

4

x

x

2

3







5

3

1

1

3

5

200

140

80

20

40

100

f x

( )

0

x

Szukamy zespolonych pierwiastków metodą Newtona - Raphsona

n

2
n

2
n

3
n

n

1

n

x

2

x

3

4

x

x

x

x













Jako punkt startowy musimy wybrać liczbę zespoloną:

np.: x

0

=i

gdzie

1

i





i

2309

.

1

8462

.

0

x

e

6056

.

3

e

1623

.

3

i

i

2

3

i

3

i

x

1

69

.

213

i

43

.

198

i

1



























x

2

=-0.50605+1.21289i

x

3

=-0.49471+1.32934i

x

4

=-0.50119+1.32409i

x

5

=-0.500059+1.322855i

x

6

=-0.5+1.32275i

błąd=-0.00083198+0.00043738i

x

d

=-0.5+1.322288i

Układy równań nieliniowych

Dany jest układ równań:













0

x

,...,

x

,

x

f

..........

..........

..........

0

x

,...,

x

,

x

f

0

x

,...,

x

,

x

f

n

2

1

n

n

2

1

2

n

2

1

1







Dla skrócenia zapisu wprowadzamy oznaczenia:





 

X

f

x

,...,

x

,

x

f

k

n

2

1

k



oraz



































)

x

,...,

x

,

x

(

f

.

.

)

x

,...,

x

,

x

(

f

)

x

,...,

x

,

x

(

f

)

X

(

F

n

2

1

n

n

2

1

2

n

2

1

1

i równanie zapisujemy krótko:

 

0

X

F



Metoda iteracji prostej

Równanie:

 

0

X

F



zapisujemy w postaci:

 

X

G

X 

i procedura iteracji prostej ma postać:





1

n

n

X

G

X





Stosowana szczególnie w przypadkach jeżeli mamy dobre

przybliżenie początkowe. Sytuacja taka występuje np. w

przypadku małej zmiany parametrów równania.

Przykład:

1

y

2

x

1

y

x

2

2

2

2









którego rozwiązaniem jest: x

1

=1; y

1

=0 oraz x

2

=-1; y

2

=0

Szukamy rozwiązania układu po małej zmianie parametrów:

1

y

01

.

2

x

95

.

0

y

x

2

2

2

2









mamy schemat iteracyjny:

01

.

2

x

1

y

y

95

.

0

x

2

1

n

n

2

1

n

n













Jako startowy punkt wybieramy: x

0

=1; y

0

=0 i mamy:

n

0

1

2

3

4

x

n

1

0.9746

8

0.97468

0.9618

34

0.96183

4

y

n

0

0

0.15771

8

0.1577

18

0.19300

6

n

5

6

7

8

x

n

0.9553

79

0.95537

9

0.95215

09

0.952151

y

n

0.1930

06

0.20834

7

0.20834

7

0.215573

6

n

9

10

11

12

x

n

0.95054

09

0.950541 0.949738

9

0.949739

y

n

0.21557

34

0.2190799

4

0.219079

7

0.220803

6

Z przedstawionych obliczeń widać, że metoda jest wolno

zbieżna i dlatego stosowana tylko w przypadkach, gdy

znamy bardzo dobrze zerowe przybliżenie. Zastosowanie

w równaniach różniczkowych.

Metoda Newtona - Raphsona

Rozwijamy funkcję f

k

(X) w szereg Taylora w otoczeniu

punktu X

i

:

0

)

x

x

(

x

f

...

)

x

x

(

x

f

)

X

(

f

...

..........

..........

..........

..........

..........

..........

..........

..........

0

)

x

x

(

x

f

...

)

x

x

(

x

f

)

X

(

f

.........

..........

..........

..........

..........

..........

..........

..........

0

)

x

x

(

x

f

...

)

x

x

(

x

f

)

X

(

f

i

,

n

n

X

X

n

n

i

,

1

1

X

X

1

n

i

n

i

,

n

n

X

X

n

k

i

,

1

1

X

X

1

k

i

k

i

,

n

n

X

X

n

1

i

,

1

1

X

X

1

1

i

1

i

i

i

i

i

i









































































Dla uproszczenia zapisu wprowadzamy macierz Jacobiego
zdefiniowaną następująco:

i

X

X

n

n

2

n

1

n

n

2

2

2

1

2

n

1

2

1

1

1

i

x

f

.

.

x

f

x

f

.

.

.

.

.

.

.

.

.

.

x

f

.

.

x

f

x

f

x

f

.

.

x

f

x

f

)

X

(

J





















































































i w postaci macierzowej możemy krótko zapisać układ równań:

   

0

X

X

J

X

F

i

i

i







gdzie oznaczono:

i

1

i

i

X

X

X









i rozwiązując symbolicznie mamy:

   

i

i

1

i

1

i

X

F

X

J

X

X









Przykład

F x y







f1 x y







f2 x y





















i równanie: F(x,y)=0

f1 x y







expx y







sin5 x



 



cosy

 





f2 x y







cosy

 

sinx

 



ln x y











Obliczamy pochodne:

p1xx y







expx y







sin5 x



 



5 expx y









cos5 x



 







p1yx y







expx y







sin5 x



 



cosy

 





p2xx y







cosx

 

cosy

 



1

x y







|

J x y







p1xx y







p2xx y







p1yx y







p2yx y





















Jakobian układu równań jest:

p2yx y







sinx

 



siny

 



1

x y







x

0

0.1



Przyjmujemy zerowe przybliżenie:

y

0

0



i liczymy z równania:

x

n 1



y

n 1















x

n

y

n













J x

n

y

n











1



F x

n

y

n













dla n=0,1,...,N

x

0

0
1
2
3
4
5
6
7
8
9

10
11

0.1

-0.141899943936
-0.029185837303
-0.044892694856
-0.036465491584
-0.035920867262
-0.035912407482
-0.035912710982
-0.035912699597
-0.035912700025
-0.035912700009
-0.035912700009



y

0

0
1
2
3
4
5
6
7
8
9

10
11

0

0.486244256751
0.743862659111
1.019360058647
1.053479318574
1.053828713899
1.053816932667
1.053817374634
1.053817358052
1.053817358674
1.053817358651
1.053817358652



Błąd:



n

F x

n

y

n











1

0.03514978332



1.191145539803



















n

F x

n

y

n











n

F x

n

y

n











10

8.754941216438

10

12





0

















n

F x

n

y

n











5

1.226398353761



10

4





5.476867793418



10

7























Obliczenia pierwiastka z pochodną liczoną numerycznie

p1xx y



h







f1 x h



y







f1 x y









h



p1yx y



h







f1 x y h









f1 x y









h



p2xx y



h







f2 x h



y







f2 x y









h



p2yx y



h







f2 x y h









f2 x y









h



Przybliżona macierz Jacobiego:

J x y



h







p1xx y



h







p2xx y



h







p1yx y



h







p2yx y



h





















Przyjmując dla kroku h:h

0.01



mamy wyniki:

x

0

0
1
2
3
4
5
6
7
8
9

10
11

0.1

-0.245620510404

0.129212408651

-0.139883972512
-0.026332867153
-0.035690080158
-0.035909711413

-0.03591266157

-0.035912699505
-0.035912700003
-0.035912700009
-0.035912700009



y

0

0
1
2
3
4
5
6
7
8
9

10
11

0

0.612829446354
0.575992105371
1.206911594234
1.055296091227
1.053366833955
1.053814123293
1.053817303719
1.053817357993
1.053817358643
1.053817358651
1.053817358652



Błąd:



1

0.541728756272



1.200733537706



















n

F x

n

y

n











1

0.03514978332



1.191145539803

















Jakobian dokładnie

Jakobian numerycznie



n

F x

n

y

n











5

1.226398353761



10

4





5.476867793418



10

7

























5

3.53469633664

10

3





1.279336142867



10

4

























n

F x

n

y

n











10

8.754941216438

10

12





0

















10

1.281363903871

10

12





1.06512021425

10

14























h

0.0001



Zmiana kroku różniczkowania

x

0

0
1
2
3
4
5
6
7
8
9

10
11

0.1

-0.243293030521

0.133173725428

-0.147314343999
-0.019012944127
-0.035661321744
-0.035912657661
-0.035912700004
-0.035912700009
-0.035912700009
-0.035912700009
-0.035912700009



y

0

0
1
2
3
4
5
6
7
8
9

10
11

0

0.597853320049
0.554664046616
1.208957700299
1.048065432775

1.05347013074

1.053817373468
1.053817358641
1.053817358652
1.053817358652
1.053817358652
1.053817358652



h=0.01

h=0.0001



1

0.510450698944



1.235991747404



















1

0.541728756272



1.200733537706



















5

3.771637969603

10

3





1.923653534576

10

5

























5

3.53469633664

10

3





1.279336142867



10

4

























10

0
0

















10

1.281363903871

10

12





1.06512021425

10

14























Metody całkowania numerycznego

Obliczana jest całka:





b

a

dx

)

x

(

f

I

Wzory Newtona - Cotesa na węzłach równoodległych

Interpolujemy funkcję f(x) za pomocą wzoru Lagrange’a:









N

0

k

k

k

N

)

x

(

f

)

x

(

)

x

(

L

gdzie

)

x

x

)...(

x

x

)(

x

x

)...(

x

x

)(

x

x

(

)

x

x

)...(

x

x

)(

x

x

)...(

x

x

)(

x

x

(

)

x

(

N

k

1

k

k

1

k

k

1

k

0

k

N

1

k

1

k

1

0

k

































Dzieląc odcinek [a,b] na N jednakowych części h:

N

a

b

h





x

a=x

0

b=x

N

h

x

k

=a+kh

Biorąc pod uwagę, że x

k

=a+kh mamy:

h

)

N

k

)...(

h

)(

h

...(

h

)

1

k

(

kh

)

Nh

a

x

)...(

h

)

1

k

(

a

x

)(

h

)

1

k

(

a

x

)...(

a

x

(

)

x

(

k





























przyjmując:

th

a

x





gdzie t[0,N] mamy:

)!

k

N

(

!

k

)

1

(

)

N

t

)...(

1

k

t

)(

1

k

t

)...(

1

t

(

t

h

)!

k

N

(

!

k

)

1

(

h

)

N

t

...(

h

)

1

k

t

(

h

)

1

k

t

...(

h

)

1

t

(

th

)

t

(

)

th

a

(

k

N

N

k

N

k

k



















































Podstawiając zamiast f(x) wielomian interpolacyjny do





b

a

dx

)

x

(

f

I

otrzymujemy:























N

0

k

k

N

0

k

k

k

N

0

k

N

0

k

k

N

dt

)

t

(

h

A

gdzie

,

)

x

(

f

A

dt

)

t

(

h

)

x

(

f

I

Błąd e

N

z jakim obliczana jest całka podaje oszacowanie:





















b

a

N

1

0

r

)

r

(

r

dx

)

x

x

)...(

x

x

)(

x

x

(

)!

1

N

(

1

C

a

,

2

)

2

/

N

(

round

2

r

gdzie

)

(

f

C

)

f

(

e

Metoda trapezów

a; 0

x; t

b; 1

f(a)

f(b)

2

h

tdt

1

1

1

h

A

2

h

dt

)

1

t

(

1

1

1

h

A

1

0

1

1

0

0





















t-1

t





)

b

(

f

)

a

(

f

2

a

b

I











)

b

(

f

)

a

(

f

2

a

b

I







błąd e

1

wynosi:

)

(

f

12

h

e

)

2

(

3

1





Metoda Simpsona

a=x

0

x

y

b=x

2





b

a

2

1

x

1





f

0

f

1

f

2

1



0



1



2

3

h

dt

)

1

t

(

t

)!

2

2

(

!

2

)

1

(

h

A

h

3

4

dt

)

2

t

(

t

)!

1

2

(

!

1

)

1

(

h

A

3

h

dt

)

2

t

)(

1

t

(

)!

0

2

(

!

0

)

1

(

h

A

2

0

2

2

2

2

0

1

2

1

2

0

0

2

0













































Całka





b

a

dx

)

x

(

f

I

obliczana metodą Simpsona jest:

 

























 







)

b

(

f

2

b

a

f

4

a

f

6

a

b

I

Błąd obliczania całki metodą Simpsona jest:





)

(

f

90

2

a

b

e

)

4

(

5

5

2











Wielomian 3-go stopnia jest całkowany dokładnie!

Kwadratury złożone Newtona-Cotesa.

x

y

a

b





b

a

2

1

x

1





Metoda trapezów

























N

1

N

1

k

k

0

N

f

f

2

f

2

h

I

Całkę





b

a

dx

)

x

(

f

I

przy podziale odcinka [a,b] na N części liczymy ze wzoru:

Obliczenia prowadzimy w schemacie z połowieniem kroku
co pozwala wykorzystać poprzednie obliczenia f

k

:











N

1

k

1

k

2

N

2

N

N

2

f

h

I

2

1

I

a

b

1
2
3
4

Błąd e

1N

złożonego wzoru trapezów przy podziale przedziału

całkowania [a,b] na N części jest:

Przykład:

4

dx

x

1

I

1

0

2











N

1

2

4

8

I

p

0.5

0.683013

0.748927

3

0.772455

I

0.785398

2

0.785398

2

0.785398

2

0.785398

2

eps

36

13

4.64

1.65

%

100

I

I

I

eps

p







  









f

N

12

a

b

e

2

3

N

1

 







 











 1

x

5

.

1

2

x

1

1

x

f

Metoda Simpsona liczba punktów podziału jest N=2

M

i całka I

N

określona jest wzorem:

































N

1

N

1

k

k

2

N

1

k

1

k

2

0

N

N

f

f

2

f

4

f

3

h

I

2
2

2

2

3

Ocena błędu:





 

 







4

4

5

N

2

f

N

180

a

b

e

4

dx

x

1

I

1

0

2











Przykład:

N

2

4

8

I

p

0.74402 0.770899

0.7802972

93

I

0.785398

2

0.785398

2

0.7853982

eps

5.3

1.85

0.65

%

100

I

I

I

eps

p





Formalna ocena błędu





 

 







4

4

5

N

2

f

N

180

a

b

e

jeszcze gorsza ze względu na czwartą pochodną.

Metoda Romberga

Ogólna metoda przyśpieszania zbieżności

Ogólnie wzór kwadratur dla całki

ma postać:

 





b

a

dx

)

x

(

f

f

I

)

x

(

f

A

)

f

(

Q

)

n

(
i

n

0

i

)

n

(
i

n







Załóżmy, że resztę można przedstawić w postaci:

)

n

c

(

O

n

c

....

n

c

n

c

)

f

(

Q

)

f

(

I

1

m

m

2

1

a

1

m

a

m

a

2

a

1

n

















gdzie

.

1

m

m

2

1

a

a

....

a

a











Niech n=sp wtedy :

)

)

sp

(

c

(

O

)

sp

(

c

....

)

sp

(

c

)

sp

(

c

)

f

(

Q

)

f

(

I

1

m

m

2

1

a

1

m

a

m

a

2

a

1

sp

















Mnożąc przez

1

a

s

i przyjmując p=n po odjęciu stronami

otrzymuje się:

)

n

c

(

O

n

c

....

n

c

n

c

)

f

(

Q

)

f

(

I

1

m

m

2

1

a

1

m

a

m

a

2

a

1

n

















)

n

c

(

O

n

c

....

n

c

n

c

)

f

(

Q

s

)

f

(

I

s

1

m

m

2

1

1

1

a

1

m

a

m

a

2

a

1

sn

a

a























)

n

c

(

O

n

c

....

n

c

1

s

)

f

(

Q

)

f

(

Q

s

)

f

(

I

1

m

m

2

1

1

a

1

1

m

a

1

m

a

1

2

a

n

sn

a



















gdzie

i

a

a

a

1

i

c

1

s

1

s

c

1

i

1









Definiując:

1

s

)

f

(

Q

)

f

(

Q

s

)

f

(

Q

1

1

a

n

sn

a

1

n







widzimy, że przybliżona wartość całki

 





b

a

dx

)

x

(

f

f

I

 

f

Q

1

n

jest obliczona z dokładnością większą bo wynoszącą











2

a

1

2

n

c

O

Ogólnie po k+1 krokach mamy:

1

s

)

f

(

Q

)

f

(

Q

s

)

f

(

Q

1

k

1

k

a

k
n

k
sn

a

1

k
n













Dla wzoru trapezów zachodzi oszacowanie:

)

n

1

(

O

n

c

...

n

c

n

c

)

f

(

Q

dx

)

x

(

f

2

m

2

m

2

m

4

2

2

1

n

b

a

















Zakładając s=2 otrzymuje się:

1

4

)

f

(

Q

)

f

(

Q

4

)

f

(

Q

1

4

)

f

(

Q

)

f

(

Q

4

)

f

(

Q

1

k

k
n

k

n

2

1

k

1

k
n

n

n

2

1

n



















Organizację obliczeń można zapisać w formie tablicy:

4

1

3
2

2
4

1

8

16

3

1

2
2

1

4

8

2

1

1

2

4

1
1

2

1

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

gdzie elementy macierzy obliczane zgodnie z podanym

powyżej algorytmem.

Uwaga - przy obliczaniu

m

2

Q korzystamy z

m

Q













m

1

i

m

2

m

2

m

m

2

)

h

)

1

i

2

(

a

(

f

h

5

.

0

Q

*

5

.

0

Q

Przykład

4

dx

x

1

I

1

0

2















5

.

0

1

1

0

1

2

1

Q

1











Q

1

=0.5

Q

2

=0.683013

1

4

Q

Q

4

Q

1

2

1
1







Q

1

1

=0.744017333

Q

4

=0.7489273

1

4

Q

Q

4

Q

2

4

1

2







Q

1

2

=0.770898733

1

4

Q

Q

4

Q

2

1
1

1

2

2

2

1







Q

2

1

=0.772690827

Q

8

=0.772455 =0.780297567

1

4

Q

1

4

Q

Q

4

Q

2

1

2

1

4

2

2
2







=0.780924156

2
2

Q

1

4

Q

Q

4

Q

3

2

1

2
2

3

3

1







=0.781054843

3

1

Q

Ocena błędu:

%

100

I

I

Q

eps

3

1





eps=0.55%

trapezy:

1.65%

Simpson: 0.65%

Przykład z ograniczoną pochodną:













2

0

2

sin

8

.

0

1

d

Q

Wstępny krok .

2

h 

 Wartość dokładna Q=2.2572051728

1

2.541602

-

-

-

-

 

2

2.2847455

921

2.1991267

89

-

-

-

 

4

2.2576215

2.2485801

36

2.251877

026

-

-

 

8

2.2572054

6

2.2570667

8

2.257632

556

2.251785

668

-

 

1
6

2.2572053

268

2.2572052

81

2.257214

514

2.257207

878

2.2572291

42

 

3
2

2.2572053

268

2.2572053

268

2.257205

329

2.257205

183

2.2572051

73

2.2572051

49

Błąd względny wynosi: 0.1E-5%.

Metoda Gaussa

Jak dobrać punkty podziału [a,b] przedziału na N części

aby uzyskać dokładność 2(N+1)?





b

a

dx

)

x

(

f

)

x

(

p

I

Rozwiązanie problemu można znaleźć za pomocą

wielomianów ortogonalnych.

Ciąg wielomianów:

)

x

(

P

),...,

x

(

P

),

x

(

P

N

1

0

jest ortogonalny w przedziale [a,b] z wagą p(x) jeżeli zachodzi:















k

=

i

dla

A

k

i

dla

0

dx

)

x

(

P

)

x

(

P

)

x

(

p

))

x

(

P

),

x

(

P

(

k

k

i

b

a

)

x

(

p

k

i

Waga p(x)=1, przedział [-1,1].

Wielomiany Legendre’a:

n

2

n

n

n

n

)

1

x

(

dx

d

!

n

2

1

)

x

(

P





 

 
 





i.t.d.

1

x

3

4

1

x

P

x

x

P

1

x

P

2

2

1

0









1 0.750.50.250 0.250.50.75 1

1

0.75

0.5

0.25

0

0.25

0.5

0.75

1

Leg0 x



(

)

Leg1 x



(

)

Leg2 x



(

)

Leg3 x



(

)

Leg4 x



(

)

x

Interpolując funkcję podcałkową f(x) za pomocą

wielomianu Legendre’a otrzymujemy:











N

0

k

k

k

1

1

)

x

(

f

A

dx

)

x

(

f

gdzie

)

x

(

P

)

x

(

P

)

2

N

(

2

A

k

'

1

N

k

2

N

k











a x

k

jest k-tym pierwiastkiem równania:

0

)

x

(

P

1

N





i błąd metody jest:

]

1

,

1

[

);

(

f

)

)!

2

N

2

)((

3

N

2

(

)

)!

1

N

((

2

e

)

2

N

2

(

3

4

3

N

2

N





















Pierwiastki i współczynniki kwadratury Gaussa z wagą p(x)=1:

N

k

1

0

0

2

2

0
1

-0.577350269189626

0.577350269189626

1
1

3

0
1
2

-0.774596669241483

0

0.774596669241483

5/9
8/9
5/9

4

0
1
2
3

-0.861136311594053
-0.339981043584856

0.339981043584856
0.861136311594053

0.347854845137454
0.652145154862546
0.652145154862546
0.347854845137454

5

0
1
2
3

4 

-0.906179845938664
-0.538469310105683

0

0.538469310105683
0.906179845938664

0.236926885056189
0.478628670499366
0.568888888888889
0.478628670499366
0.236926885056189

węzły x

k

współczynniki A

k

W przypadku całkowania w przedziale [a,b] dokonujemy
jego transformacji na przedział [-1,1]:

  



1

,

1

b

,

a









dt

2

a

b

dx

1

t

2

a

b

a

x













Przykład

4

dx

x

1

I

1

0

2











n

1

2

3

4

5

I

n

0.86602

54

0.79611

301

0.78725

969

0.78705

74

0.78629

544

błąd
%

10.3

1.36

0.237

0.211

0.114

Ta sama całka, ale liczona w ten sposób, że przedział całkowania

jest podzielony na 5 równych części i w każdym przedziale

zastosowano aproksymację N=1

I

1

=0.792996956 z błędem 0.97%

Całki typu:



 

  











b

a

dx

x

f

a

x

x

b

Transformacja przedziału całkowania:

  



1

,

1

b

,

a









dt

2

a

b

dx

1

t

2

a

b

a

x













Pozwala zapisać całkę w postaci:



 

  













1

1

0

dt

t

f

t

1

t

1

A

gdzie



















 



1

0

2

a

b

A

Wystarczy więc rozważyć całkę:



 

  















1

1

dt

t

f

t

1

t

1

I

Mamy obliczyć całkę funkcji f(t) z wagą

Wielomiany ortogonalne na przedziale [-1,1] z wagą w(t)

są nazywane wielomianami Jacobiego i mają postać:





    

 





 







t

1

t

1

dt

d

t

1

t

1

!

n

2

1

t

P

n

n

n

n

n

n

,

n

































i całkę I można obliczyć z zależności:



 

  

 



















N

0

k

k

k

1

1

t

f

A

dt

t

f

t

1

t

1

  

 













t

1

t

1

t

w

gdzie



 

 





 

 







 





 

k

,

2

N

k

,

1

N

k

t

P

t

P

2

N

2

N

!

2

N

2

N

2

N

4

N

2

2

A









































































a t

k

jest pierwiastkiem równania:





 

0

t

P

k

,

1

N









Przypadki szczególne:

i

 

 













0

1

x

dt

t

exp

t

x





!

n

1

n







dla n=0,1,2,...





























 



n

2

1

n

2

5

3

1

2

1

n



dla n=1,2,3,...

 





 n

dla n=0,-1,-2,...

Typowe funkcje wagowe

.

1 ==-0.5 czyli

 

2

t

1

1

t

w





a całka ma postać:

 









1

1

2

dt

t

1

t

f

I

W tym przypadku wielomiany Jacobiego są powiązane

z wielomianami Czebyszewa T

n

(t) związkiem:





 

 

t

T

C

t

P

n

n

5

.

0

,

5

.

0

n







gdzie C

n

– stała , a

 

 





t

arccos

n

cos

t

T

n





a wykres:

  

 













t

1

t

1

t

w

1 0.750.50.250 0.250.50.75 1

1.25

0.94

0.63

0.31

0

0.31

0.63

0.94

1.25

Jac 0 0.5





0.5





x



(

)

Jac 1 0.5





0.5





x



(

)

Jac 2 0.5





0.5





x



(

)

Jac 3 0.5





0.5





x



(

)

Jac 4 0.5





0.5





x



(

)

x

Łatwo więc znaleźć pierwiastki t

k

z równania:

 





0

t

arccos

N

cos

t

T

k

k

N







i mamy:



















N

2

1

k

2

cos

t

k

dla k=1,2,3,...,N

Korzystając z odpowiednich tożsamości dla funkcji Czebyszewa
i Jacobiego można wykazać, że współczynniki wagi są:





 

2

2

N

2

N

k

N

C

N

!

N

5

.

0

N

2

A









a ponieważ nie są zależne od k, więc dla danego N powinny być
stałe i równe A. Stałą A łatwo wyznaczymy jeżeli przyjmiemy, że
f(t)=1 i mamy:

NA

A

dt

t

1

N

1

k

k

1

1

2



















i stąd

N

A 



mamy ostatecznie:

 

 

f

R

N

2

1

k

2

cos

f

N

dt

t

1

t

f

N

N

1

k

1

1

2









































gdzie R

N

(f) jest błędem aproksymacji wynoszącym:

 

 

 

 

!

N

2

f

2

f

R

N

2

1

N

2

N











i –1<<1

Przykład













0

2

2

cos

k

1

d

I

gdzie |k|<1

Podstawiamy: cos=t i

2

t

1

dt

d







czyli:

 











1

1

2

2

t

1

kt

1

dt

I

w tym przypadku:

 











1

1

2

2

t

1

kt

1

dt

I

funkcja f(t) ma postać:

 

 

2

kt

1

1

t

f





i na mocy zależności:

 

 

f

R

N

2

1

k

2

cos

f

N

dt

t

1

t

f

N

N

1

k

1

1

2









































mamy:







































N

1

i

2

2

0

2

2

N

2

1

i

2

cos

k

1

1

N

cos

k

1

d

I

0

0.2

0.4

0.6

0.8

1

2

4

6

8

I k

( )

Ip k 2



(

)

Ip k 4



(

)

Ip k 6



(

)

k

0

0.2

0.4

0.6

0.8

1

20

6.67

33.33

60

eps k 2



(

)

eps k 4



(

)

eps k 6



(

)

eps k 8



(

)

k

%

100

)

k

(

I

)

n

,

k

(

Ip

)

k

(

I

)

n

,

k

(

eps





2. ==0.5

Funkcja wagowa:

 

2

t

1

t

w





 









1

1

2

dt

t

f

t

1

I

a całka ma postać:

W tym przypadku funkcje Jacobiego można wyrazić przez
funkcje Czebyszewa II – go rodzaju U

n

(t):





 









 

t

U

!

1

n

!

n

2

!

1

n

2

t

P

n

n

2

5

.

0

,

5

.

0

n







gdzie

 





 





2

n

t

1

t

arccos

1

n

sin

t

U







wykres wielomianu Jacobiego dla ==0.5 jest:

1 0.750.50.250 0.250.50.75 1

2.5

1.88

1.25

0.63

0

0.63

1.25

1.88

2.5

Jac 0 0.5



0.5



x



(

)

Jac 1 0.5



0.5



x



(

)

Jac 2 0.5



0.5



x



(

)

Jac 3 0.5



0.5



x



(

)

Jac 4 0.5



0.5



x



(

)

x

 

 

f

R

1

N

k

cos

f

1

N

k

sin

1

N

dt

t

f

t

1

N

N

1

k

2

1

1

2



























































gdzie błąd określa wzór:

 

 

 

 

!

N

2

f

2

f

R

N

2

N

2

N







-1<<1

Przykład:















1

0

2

2

dx

4

x

4

x

x

x

2

I

transformacja przedziału [0,1] do przedziału [-1,1]:1

t

x





i całka przyjmuje postać:

















1

1

2

2

dt

1

t

2

t

t

1

I

Z oceny błędu widać, że wystarczy wybrać N=3, a wynik
będzie dokładny i mamy: I=1.96349541

a wzór dla całki jest:

Jeszcze jeden przypadek szczególny

3. =0.5 =-0.5

Funkcja wagowa jest:

 

t

1

t

1

t

w







i liczymy całkę:

 











1

1

dt

t

f

t

1

t

1

I

Również w tym przypadku można wyrazić wielomiany
Jacobiego przez wielomiany Czebyszewa:





 

 

 

 

 

1

t

t

T

t

T

!

n

2

!

n

2

t

P

n

1

n

2

n

2

5

.

0

,

5

.

0

n











i wykres jest:

1 0.750.50.250 0.250.50.75 1

2.5

1.88

1.25

0.63

0

0.63

1.25

1.88

2.5

Jac 0 0.5



0.5





x



(

)

Jac 1 0.5



0.5





x



(

)

Jac 2 0.5



0.5





x



(

)

Jac 3 0.5



0.5





x



(

)

Jac 4 0.5



0.5





x



(

)

x

i kwadraturę można zapisać:

 

 

f

R

1

N

2

k

2

cos

f

1

N

2

k

sin

1

N

2

4

dt

t

f

t

1

t

1

N

N

1

k

2

1

1





























































gdzie

 

 

 

 







N

2

N

2

N

f

!

N

2

2

f

R

gdzie -1<<1

Całki na przedziale nieograniczonym

Można wybrać specjalne reguły całkowania za pomocą

wielomianów ortogonalnych Laguerre’a lub Hermite’a

w przedziale [0,] lub [-,+] odpowiednio.

 

 





















dx

x

f

I

dx

x

f

I

2

a

1

Całki postaci:

 









0

x

dx

x

f

e

x

gdzie >-1

wagą jest

 

x

e

x

x

w







Wielomianami ortogonalnymi na półosi [0,) z wagą w(x)

są wielomiany Laguerre’a:

 

   





x

n

n

n

x

n

n

e

x

dx

d

e

x

1

x

L

















Dla =0 mamy:

 

 













N

1

k

k

k

0

x

x

f

A

dx

x

f

e

a x

k

– zera wielomianu Laguerre’a i współczynniki A

k

podaje

tablica:

N=1 x

1

=1.00000000

A

1

=1.00000000

N=2 x

1

=0.585786437

627
x

2

=3.414213562

373

A

1

=0.853553390593

A

2

=0.146446609407

N=3 x

1

=0.215774556

783
x

2

=2.294280360

279
x

3

=6.289945082

937

A

1

=0.711093009929

A

2

=0.278517733569

A

3

=0.0103892565016

N=4 x

1

=0.322547689

619
x

2

=1.745761101

158
x

3

=4.536620296

921
x

4

=9.395070912

301

A

1

=0.603154104342

A

2

=0.357418692438

A

3

=0.038887908515

A

4

=0.539294705561·1

0

-3

N=
5

x

1

=0.2635603197

18
x

2

=1.4134030591

07
x

3

=3.5964257710

41
x

4

=7.0858100058

59
x

5

=12.640800844

276

A

1

=0.521755610583

A

2

=0.398666811083

A

3

=0.0759424496817

A

4

=0.00361175867992

A

5

=233699723858·10

-4

N=

6

x

1

=0.2228466041

79
x

2

=1.1889321016

73
x

3

=2.9927363260

59
x

4

=5.7751435691

05
x

5

=9.8374674183

83
x

6

=15.982873980

602

A

1

=0.458964673950

A

2

=0.417000830772

A

3

=0.113373382074

A

4

=0.0103991974531

A

5

=0.261017202815·1

0

-3

A

6

=0.89854790643·10

-6

Przykład:
Dana jest całka:













0

x

2

x

dx

e

1

xe

I

Dokładna wartość całki jest:

8

I

2





Obliczenia metodą trapezów przy obcięciu górnej granicy:

 











x

0

t

2

t

dt

e

1

te

x

In

Błąd obliczeń liczymy z zależności:

 

100

I

I

x

In

)

x

(

err





Wyniki przedstawia tabela:

x=1

-61.6%

x=5

-3.28%

x=10

-0.04%

x=100

-2.28·10

-

10

%

x=1000

2.26·10

-9

%

dlaczego?

Obliczenia korzystając z wielomianu Laguerre’a

N=1

-6.26%

N=2

-0.7%

N=3

0.068%

N=4

0.048%

N=5

7.82·10

-3

%

N=6

-2.52·10

-3

%

Można też zastosować następujące postępowanie dla całki:

 

0

a

;

dx

x

f

I

a

1











podstawiamy:

t

dt

dx

e

t

x









i mamy:

 









a

e

0

t

1

t

dt

e

f

I

Warunkiem istnienia całki I

1

jest:

 

0

t

e

f

t

0

t

lim

















a więc stosując np. metodę trapezów trzeba w punkcie t=0

uwzględnić wartość graniczną.

Całkę

 













dx

x

f

I

2

najwygodniej rozbić na dwie:

 

 

 





























0

0

2

dx

x

f

dx

x

f

dx

x

f

I

które liczymy jak poprzednio.

Funkcja podcałkowa posiada osobliwość wewnątrz

przedziału całkowania

Np.:







2

0

x

1

dx

I

ogólnie:

 





b

a

dx

x

f

I

i w punkcie x

0

[a,b] występuje osobliwość.

1. Jeżeli istnieje

   

0

x

x

x

f

x

f

lim

0





to wydzielamy osobliwość:

 

   







  

0

b

a

0

b

a

x

f

b

a

dx

x

f

x

f

dx

x

f

I















Przykład:















2

2

2

2

dx

5

x

6

x

4

x

3

x

I

wewnątrz przedziału całkowania mianownik ma miejsce
zerowe x

0

=1, a granica funkcji w tym punkcie wynosi:

25

.

1

5

x

6

x

4

x

3

x

lim

2

2

1

x















czyli

4

25

.

1

dx

25

.

1

5

x

6

x

4

x

3

x

dx

5

x

6

x

4

x

3

x

I

2

2

2

2

2

2

2

2











































5

dx

25

.

1

5

x

6

x

4

x

3

x

I

2

2

2

2



























w obliczeniach numerycznych mamy funkcję podcałkową:

 

























1

x

dla

0

1

x

dla

25

.

1

5

x

6

x

4

x

3

x

x

f

2

2

z której całkę w przedziale [-2,2] liczymy dowolną metodą.

2. Funkcja f(x) posiada słabą osobliwość w punkcie x

0

[a,b]:

 







0

x

x

x

g

)

x

(

f

gdzie 0<<1 oraz g(x

0

)0.

i mamy:

 

   

 



























b

a

b

a

0

0

0

0

b

a

0

x

x

dx

x

g

dx

x

x

x

g

x

g

dx

x

x

x

g

I

Ostatnią całkę można obliczyć analitycznie:

 

 



























1

x

a

x

b

x

g

x

x

dx

x

g

1

0

1

0

0

b

a

0

0

Natomiast funkcja podcałkowa:

 

   









0

0

x

x

x

g

x

g

x

f

ma w punkcie x

0

granicę:

 

0

x

f

lim

0

x

x





jeżeli







0

x

x

dx

dg

Równanie różniczkowe rzędu n-go

 





0

x

,...,

x

,

x

,

t

f

n





przykład:

























t

sin

E

x

x

x

1

a

x

2

2





Układ równań różniczkowych I-go rzędu w postaci
normalnej:













n

2

1

n

n

n

2

1

2

2

n

2

1

1

1

x

,...,

x

,

x

,

t

f

x

........

..........

..........

x

,...,

x

,

x

,

t

f

x

x

,...,

x

,

x

,

t

f

x













Numeryczne metody całkowania równań różniczkowych

Sprowadzenie równania rzędu n do układu równań na
przykładzie:

























t

sin

E

x

x

x

1

a

x

2

2





oznaczamy:

 

 

2

2

2

1

1

x

x

dt

d

x

dt

d

x

x

x

x

x

x























i ostatecznie:









1

2

2

2

1

2

2

1

x

x

x

1

a

t

sin

E

x

x

x























Równanie różniczkowe rzędu n może praktycznie zawsze

być zapisane w formie układu n równań I-go rzędu

mówimy, że mamy układ równań w postaci normalnej.

Zagadnienie początkowe dla równania różniczkowego:

Przykład:

 

t

f

x

x









całka ogólna:

0

x

x













t

exp

A

x







Rozwiązanie:

Całka szczególna metodą uzmienniania stałej:



  

t

f

t

exp

A









czyli

 

 













t

0

d

exp

f

C

A

Rozwiązanie równania:

 

t

f

x

x









ma postać:





 































t

0

d

t

exp

f

t

exp

C

x

Stałą C wyznaczamy z warunków początkowych

Dla równania różniczkowego I-go rzędu potrzebujemy
jeden warunek początkowy:





0

x

0

t

x





i dla wyznaczenia stałej C mamy równanie:

C

x

0



i rozwiązanie ma postać:

 





 































t

0

0

d

t

exp

f

t

exp

x

t

x

Przykład równania różniczkowego II rzędu:

R

0

E

C

R

L

t=0

i(t)

u

c

c

u

C

i





0

u

u

RC

u

LC

c

c

c











Warunki początkowe:

 

E

R

R

R

0

u

0

c





0

R

R

E

)

0

(

i







 

E

R

R

R

0

u

0

c





0

R

R

E

)

0

(

i







Biorąc pod uwagę, że

c

u

C

i



zapisujemy warunki w postaci:

 





C

R

R

E

u

E

R

R

R

0

u

0

0

t

c

0

c















Warunków początkowych należy postawić tyle i nie więcej

ile wynosi rząd równania różniczkowego

0

u

u

RC

u

LC

c

c

c











 





C

R

R

E

u

E

R

R

R

0

u

0

0

t

c

0

c















Rozwiązanie:

Niech równanie charakterystyczne:

0

1

RC

LC

2











ma dwa różne pierwiastki rzeczywiste

2

1

,



wtedy całka ogólna ma postać:

 

t

2

t

1

c

2

1

e

A

e

A

t

u









Stałe A

1

i A

2

wyznaczamy z warunków początkowych:

 





C

R

R

E

u

E

R

R

R

0

u

0

0

t

c

0

c















czyli





C

R

R

E

A

A

E

R

R

R

A

A

0

2

2

1

1

0

2

1



















Rozwiązując powyższy układ równań wyznaczamy stałe

a następnie znajdujemy napięcie na kondensatorze.

Podsumowanie:

1. Mamy kłopoty z rozwiązywaniem równań o stałych
współczynnikach, jeżeli rząd wyższy od dwóch
równanie charakterystyczne jest wielomianem
rzędu n i nie znamy wzorów na pierwiastki.
2. Równań o zmiennych w czasie współczynnikach
i nieliniowych nie potrafimy rozwiązać w ogólnym
przypadku często już dla rzędu I-go

Konieczność użycia metod numerycznych

Ze względu na łatwość zbudowania ogólnego algorytmu

obliczenia prowadzimy dla układu równań różniczkowych

I-go rzędu w postaci normalnej

który po wprowadzeniu formalnego zapisu w postaci wektorów:

 

 

 

 

























t

x

t

x

t

x

t

X

n

2

1



 

 

 

 

























t

x

t

x

t

x

t

X

n

2

1























n

2

1

n

n

n

2

1

2

2

n

2

1

1

1

x

,...,

x

,

x

,

t

f

x

........

..........

..........

x

,...,

x

,

x

,

t

f

x

x

,...,

x

,

x

,

t

f

x





















































n

2

1

n

n

2

1

2

n

2

1

1

x

,...,

x

,

x

,

t

f

x

,...,

x

,

x

,

t

f

x

,...,

x

,

x

,

t

f

X

,

t

F



i układ zapisujemy:





X

,

t

F

X 







X

,

t

F

X 



Wniosek – Jeżeli opanujemy metodę numerycznego
rozwiązywania jednego równania różniczkowego
pierwszego rzędu, czyli

 

x

,

t

f

x 



to wyniki łatwo uogólnimy na układ n równań
I-go rzędu w postaci normalnej

Rozpoczynamy od następującego zadania:

 

x

,

t

f

x 



z warunkiem początkowym:





0

x

0

t

x





Należy znaleźć rozwiązanie dla t[0,T]

Metoda aproksymacji wielomianowej

metody wielokrokowe

 

t

,

x

f

x 



Przyjmujemy algorytm w postaci























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

t

k

=kh

Algorytm nazywamy

jawnym

, jeżeli sumowanie rozpoczyna

się od 0 w przeciwnym przypadku mówimy, że algorytm jest

niejawny

Wyznaczenie współczynników a

i

, b

i

.

Metoda jest dokładna jeżeli rozwiązaniem równania:

 

t

,

x

f

x 



jest wielomian stopnia zerowego, czyli równanie ma postać:

0

x

.



którego rozwiązaniem jest:

 

0

c

t

x









p

1

i

0

i

0

c

a

c

podstawiając do























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

mamy

ponieważ f(x,t)=0. Dzieląc przez c

0

Wyznaczenie współczynników a

i

, b

i

.

Metoda jest dokładna dla wielomianu stopnia zerowego
jeżeli rozwiązaniem równania:

 

t

,

x

f

x 



jest wielomian stopnia zerowego, czyli równanie ma postać:

0

x

.



którego rozwiązaniem jest:

 

0

c

t

x









p

1

i

0

i

0

c

a

c

podstawiając do























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

mamy

ponieważ f(x,t)=0. Dzieląc przez c

0

otrzymujemy warunek:

1

a

p

0

i

i







Jeżeli współczynniki a

i

dobierzemy tak, aby spełnić

powyższy warunek, to algorytm jest dokładny dla
wielomianów stopnia zerowego.

Żądamy, jeżeli rozwiązaniem równania

0

1

c

t

c

x





 

t

,

x

f

x 



jest wielomian pierwszego stopnia

to algorytm























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

jest dokładny.

Równanie spełniane przez wielomian

0

1

c

t

c

x





ma postać:

1

.

c

x 

czyli

c

)

t

,

x

(

f

1



a więc algorytm:























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

uwzględniając, że t

i

=ih przyjmuje postać:























p

1

i

i

i

p

0

i

0

1

i

0

1

c

b

h

]

c

h

)

i

n

(

c

[

a

c

h

)

1

n

(

c

Biorąc pod uwagę, że dla wielomianu stopnia zerowego
mamy warunek:

1

a

p

0

i

i







Z równania:























p

1

i

i

i

p

0

i

0

1

i

0

1

c

b

h

]

c

h

)

i

n

(

c

[

a

c

h

)

1

n

(

c

otrzymujemy:

1

b

ia

p

1

i

i

p

0

i

i

















Powtórzmy jeszcze jako ćwiczenie rozumowanie dla
wielomianu II-go stopnia:

0

1

2

2

c

t

c

t

c

x







który spełnia równanie

1

2

.

c

t

c

2

x





czyli

1

2

c

t

c

2

)

t

,

x

(

f





Przyjmując dla skrócenia zapisu t

n

=0 algorytm:























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

zapisujemy:

i uwzględniając:

1

2

i

n

i

n

0

1

2

2

i

n

0

1

2

2

1

n

c

)

ih

(

c

2

)

t

,

x

(

f

c

)

ih

(

c

)

ih

(

c

x

c

h

c

h

c

x

























































p

1

i

1

2

i

p

0

i

0

1

2

2

i

0

1

2

2

]

c

ih

c

2

[

b

h

]

c

ih

c

)

ih

(

c

[

a

c

h

c

h

c

a biorąc pod uwagę dwa poprzednie warunki:

1

a

p

0

i

i







1

b

ia

p

1

i

i

p

0

i

i

















otrzymujemy:













1

i

b

2

i

a

p

1

i

i

2

i

Dla wielomianu stopnia k:

0

1

1

k

1

k

k

k

c

t

c

...

t

c

t

c

x















równie różniczkowe ma postać

1

2

k

1

k

1

k

k

.

c

...

t

c

)

1

k

(

t

kc

x

















Przyjmując podobnie jak dla drugiego stopniaih

t

i



































k

1

m

1

m

m

i

n

i

n

k

0

m

m

m

i

n

k

0

m

m

m

1

n

)

ih

(

mc

)

t

,

x

(

f

)

ih

(

c

x

h

c

x

i podstawiając do:























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

otrzymujemy:

 

 

























p

1

i

k

1

m

1

m

m

i

p

0

i

k

0

m

m

m

i

k

0

m

m

m

)

ih

(

mc

b

h

)

ih

(

c

a

h

c

Porównując współczynniki przy tych samych potęgach
h i poprzednie k-1 równań znajdujemy:

1

)

i

(

b

k

)

i

(

a

p

1

i

1

k

i

p

0

i

k

i





















Algorytm wielokrokowy:























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

jest dokładny dla wielomianu stopnia k, jeżeli współczynniki
a

i

, b

i

spełniają następujący układ k równań:

1

a

p

0

i

i







1

b

ia

p

1

i

i

p

0

i

i





























1

i

b

2

i

a

p

1

i

i

2

i

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

1

)

i

(

b

k

)

i

(

a

p

1

i

1

k

i

p

0

i

k

i





















Każdy algorytm spełniający warunki dla k>1 nazywamy

zwartym. Jeżeli metoda jest dokładna dla wielomianu stopnia

k, to mówimy o metodzie rzędu k-go.

Liczba niewiadomych do wyznaczenia w algorytmie:























p

0

i

p

1

i

i

n

i

n

i

i

n

i

1

n

)

t

,

x

(

f

b

h

x

a

x

Współczynników a

i

wynosi: p+1

Współczynników b

i

wynosi: p+2

Całkowita liczba niewiadomych: 2p+3

1

)

i

(

b

k

)

i

(

a

p

1

i

1

k

i

p

0

i

k

i





















Układ równań:

k=0,1,...,n

można spełnić pod warunkiem, że liczba równań n+1  2p+3

W metodzie Adamsa - Bashfortha przyjmujemy p=n-1

Liczba niewiadomych wynosi: 2(n-1)+3=2n+1, więc przy

spełnieniu n+1 równań możemy dowolnie wybrać

n współczynników

Przyjmujemy n-1 współczynników a

1

=a

2

=....=a

n-1

=0

Pierwsze

z równań przy spełnieniu powyższych warunków

przyjmuje postać:

1

a

p

0

i

i







1

a

0



Jako n-ty dowolnie wybrany współczynnik przyjmujemy b

-1

=0

A więc schemat Adamsa - Bashfortha jest

schematem jawnym

Dla n=1 mamy p=0 i równanie wyznaczające b

0

1

)

i

(

b

k

)

i

(

a

p

1

i

1

k

i

p

0

i

k

i





















jest

1

b

0



Przy n=2 będzie p=1 i mamy dwie niewiadome b

0

i b

1

spełniające układ równań:

2

1

b

1

b

b

1

1

0









czyli

2

1

b

;

2

3

b

1

0



































































































































n

1

.

3

1

2

1

1

b

.

b

b

b

))

1

n

(

(

...

)

2

(

)

1

(

0

.

...

.

.

.

)

1

n

(

...

4

1

0

)

1

n

(

...

2

1

0

1

...

1

1

1

1

n

2

1

0

1

n

n

n

2

Ogólnie jeżeli mamy n niewiadomych współczynników b

i

,

to wyznaczamy je z równań

Rozwiązanie powyższego układu równań przedstawia

tabela:

Współczynniki b

k

metoda Adamsa -

Bashfortha

n

b

0

b

1

b

2

b

3

b

4

b

5

1

1

2

3/2

-1/2

3

23/12

-16/12

5/12

4

55/24

-59/24

37/24

-9/24

5

1901/72

0

-

2774/720

2616/72

0

-

1274/720

251/720

6

4277/14

40

-

7923/144

0

9982/14

40

-

7298/144

0

2877/14

40

-475/1440

Algorytm Adamsa - Bashfortha rzędu pierwszego jest

nazywany algorytmem Eulera i ma postać:





i

i

i

1

i

t

,

x

hf

x

x







Przykład:

 

t

30

sin

10

x

3

x







z warunkiem początkowym x(0)=0

Zapisujemy równanie w postaci normalnej:

 

x

3

t

30

sin

10

x







czyli w tym przypadku funkcja f(x,t) jest

 

 

x

3

t

30

sin

10

t

,

x

f





Wybór kroku całkowania h:

Wybór kroku całkowania h:

Równanie jednorodne ma postać:

0

x

3

x







Całka ogólna tego równania ma postać:





t

3

exp

A

x





Przebieg rozwiązania jest scharakteryzowany przez

wielkość tłumienia, które w tym przypadku wynosi

a=3 i charakterystyczny czas wynosi 1/a=1/3 s.

Krok czasowy h należy wybrać co najwyżej

h<1/(10a), co w tym przypadku pozwala wybrać

h=0.03s.

Należy również uważnie rozważyć funkcję wymuszającą,
która w rozpatrywanym przypadku ma postać:

 

t

30

sin

10

Okres T zmian analizowanej funkcji wynosi:

s

21

.

0

30

2

T







Rysunek funkcji sinus jest „gładki”, jeżeli podzielimy

okres na co najmniej 20 kroków, czyli krok h powinien

wynosić h<T/20, a więc h=0.01s.

t

0

h

2h

3h

4h

Z obu ograniczeń h=0.03s i h=0.01s wybieramy mniejszy
czyli w tym przypadku przyjmujemy h=0.01s.
Algorytm przy przyjętym kroku h ma postać:









n

n

1

n

x

3

n

01

.

0

30

sin

10

01

.

0

x

x













x

0

=0

i mamy:









02955

.

0

0

3

01

.

0

30

sin

10

01

.

0

x

01

.

0

t

1

1















następny:









08513

.

0

02955

.

0

3

02

.

30

sin

10

01

.

02955

.

0

x

02

.

0

t

2

2















itd..

0 15 30 45 60 75 90105120135150

0.5

0.25

0

0.25

0.5

0.75

1

x

n

y n h



(

)

n

h=0.01s

y(t) rozwiązanie dokładne

0 15 30 456075 90105120135150

0.1

0.067

0.033

0

0.033

0.067

0.1

x

n

y n h



(

)



n

0 3 6 9 12 15 18 21 24 27 30

0.5

0.25

0

0.25

0.5

0.75

1

x

n

y n h



(

)

n

h=0.05s

0 3 6 9 12 1518 2124 27 30

0.5

0.33

0.17

0

0.17

0.33

0.5

x

n

y n h



(

)



n

0 1.5 3 4.5 6 7.5 9 10.51213.515

1

0.67

0.33

0

0.33

0.67

1

x

n

y n h



(

)

n

h=0.1s

0 1.5 3 4.5 6 7.59 10.51213.515

1

0.67

0.33

0

0.33

0.67

1

x

n

y n h



(

)



n

Równanie drugiego rzędu na przykładzie wahadła
matematycznego o długości l:

0

sin

l

g

dt

d

2

2











Wprowadzamy nowe zmienne
celem zastąpienia układem
równań I-go rzędu:















sin

l

g

dt

d

dt

d

Przyjmujemy wahadło o długości l=1m i g=10m/s

2

Warunki początkowe:

 

 

0

dt

d

0

18

0

0

t

















Jeżeli zlinearyzujemy równanie:

0

sin

l

g

dt

d

2

2









to przyjmie on postać:

0

l

g

dt

d

2

2









którego całka ogólna ma postać:































t

l

g

cos

A

t

l

g

sin

A

2

1

Okres T drgań zlinearyzowanego wahadła wynosi:

s

2

g

l

2

T







czyli krok należy przyjąć

1

.

0

h

20

T

h





ponieważ mamy nieliniowe równanie, więc dla bezpieczeństwa
przyjmujemy h=0.001s.

Jawny schemat Eulera dla układu równań















sin

l

g

dt

d

dt

d

przyjmuje postać:

n

n

1

n

n

n

1

n

sin

l

g

h

h

























z warunkiem startowym:

0

18

0

0











0

4000 8000

1.2.10

4

1.6.10

4

2.10

4

0.2

0.13

0.067

0

0.067

0.13

0.2



n

x n h



(

)

n

h=0.001

0

4000 80001.2.10

4

1.6.10

4

2.10

4

0.05

0.03

0.01

0.01

0.03

0.05



n

x n h



(

)



n

błąd

0.20.12

0.04

0.04

0.120.2

1

0.6

0.2

0.2

0.6

1

y1 

 

y2 

 



ω

Płaszczyzna fazowa – rozwiązanie analityczne:

 





0

cos

cos

l

g

2

dt

d



















0.2 0.12 0.04 0.04 0.12 0.2

1

0.6

0.2

0.2

0.6

1



n



n

płaszczyzna fazowa

0

400 800 1200 1600 2000

0.5

0.33

0.17

0

0.17

0.33

0.5



n

x n h



(

)

n

h=0.01

0

400 800 1200 1600 2000

0.5

0.3

0.1

0.1

0.3

0.5



n

x n h



(

)



n

błąd

0.5

0.3

0.1 0.1

0.3

0.5

2

1.2

0.4

0.4

1.2

2



n



n

płaszczyzna fazowa

0

200 400 600 800 1000

200

100

0

100

200

300

400



n

x n h



(

)

n

h=0.05

0

200 400 600 800 1000

200

80

40

160

280

400



n

x n h



(

)



n

błąd

200 80

40

160 280 400

10

4

2

8

14

20



n



n

płaszczyzna fazowa