Index of /rozprawy2/10512

Pełen tekst

(1)AGH Akademia Górniczo - Hutnicza im. St. Staszica w Krakowie Wydział Odlewnictwa Katedra Inżynierii Procesów Odlewniczych. Rozprawa doktorska Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. Paweł Leszek Żak. Promotor: Prof. dr hab. inż. Józef Szczepan Suchy. Kraków 2012.

(2) Składam serdeczne podziękowania Władzom Wydziału Odlewnictwa Akademii Górniczo-Hutniczej w Krakowie za przeprowadzenie przewodu doktorskiego. Dziękuję Komisji Europejskiej za wsparcie finansowe w latach 2007 - 2010, udzielone w ramach projektu Marie Curie TOK-DEV MTKD-CT-2006-042468 pt. Development of environmentally friendly cast alloys and composites - CastModel. Ministerstwu Nauki i Szkolnictwa Wyższego za wsparcie finansowe udzielone w ramach projektów: Projekt badawczy: N507-44-66-34 pt. Modelowanie mikro makro krzepnięcia układu heterofazowego na przykładzie kompozytu na osnowie stopu AZ91. Grant Dziekański: 15.11.170.424 pt. Opracowanie algorytmu rozwiązywania nieliniowego równania różniczkowego transportu ciepła przy pomocy metody Kwadratur Różniczkowych Sterowanego Rzędu.. Dziękuję pracownikom Wydziału Odlewnictwa AGH, którzy przyczynili się do powstania tej pracy, a w szczególności promotorowi prof. dr hab. inż. Józefowi Szczepanowi Suchemu oraz dr inż. Januszowi Lelito.

(3) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. SPIS TREŚCI: Rozdział 1: Podstawy teoretyczne kwadratury różniczkowej oraz matematycznego opisu procesów przewodzenia ciepła 5 1.1. Podstawy matematyczne opisu przewodzenia ciepła. 7. 1.1.1. Podstawowe równania rządzące procesami migracji ciepła. 7. 1.1.2. Opis warunków jednoznaczności. 7. 1.1.2.1. Warunki geometryczne. 8. 1.1.2.2. Warunki fizyczne. 8. 1.1.2.3. Warunki początkowe. 8. 1.1.2.4. Warunki brzegowe. 8. 1.1.3. Równanie różniczkowe Fouriera - Kirchhoffa dla obszarów jednowymiarowych 9. 1.1.3.1. Przypadek szczególny: równanie transportu ciepła bez funkcji źródła. 9. 1.1.4. Funkcja udziału objętości zakrzepłej. 10. 1.1.5. Zastępcza pojemność cieplna strefy dwufazowej. 11. 1.1.6. Podsumowanie. 12. 1.2. Opis metody kwadratur różniczkowych. 12. Definicja kwadratury różniczkowej. 13. 1.2.1. 1.2.2 Wyznaczanie współczynników wagowych metody KR pierwszego rzędu. dla. pochodnej 14. 1.2.2.1. Baza kanoniczna przestrzeni wielomianów. 14. 1.2.2.2. Jednoznaczność współczynników wagowych KR dla różnych baz [1]. 15. 1.2.2.3. Baza wielomianów interpolacyjnych Lagrange’a. 16. 1.2.3 Wyznaczanie współczynników metody KR dla pochodnych drugiego i wyższych rzędów 18 1.2.3.1. Metoda oparta o własności macierzy współczynników wagowych. 18. 1.2.3.2 Metoda rekurencyjnego wyznaczania współczynników wagowych dla pochodnych drugiego i wyższych rzędów 18 1.3. Dokładność przybliżenia metodą KR. 1.3.1. Zagadnienie najlepszej aproksymacji. 20 20. 1.3.1.1. Zagadnienie najlepszej aproksymacji w przestrzeni wektorowej. 20. 1.3.1.2. Typy aproksymacji. 21 Strona 1.

(4) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.3.2. Interpolacja wielomianowa. 1.3.2.1. 22. Dokładność interpolacji. 23. 1.3.3 Teoretyczna analiza dokładności przybliżenia pochodnej metodą kwadratur różniczkowych 27. 1.4. 1.3.3.1. Przybliżenie funkcji wielomianem interpolacyjnym. 28. 1.3.3.2. Związek KR z pochodną wielomianu interpolacyjnego. 29. 1.3.3.3. Przybliżenie dowolnej pochodnej funkcji kwadraturą różniczkową. 29. Udoskonalone metody KR. 31. 1.4.1 Równoważność metody kwadratur różniczkowych i schematu różnicowego najwyższego rzędu 31 1.4.2 Metoda kwadratur różniczkowych zmiennego rzędu (ang. Variable Order Differential Quadrature) 33 1.4.3 Metoda kwadratur różniczkowych Adaptative Differential Quadrature). ograniczonego. zasięgu. (ang. Local 34. 1.4.4 Metoda umiejscowionej kwadratury różniczkowej (ang. Localized Differential Quadrature Method) 35 1.5. Podsumowanie. Rozdział 2:. 36. Cel i tezy pracy. 38. 2.1. Cel pracy. 38. 2.2. Tezy pracy. 38. Rozdział 3: Metoda kwadratur różniczkowych sterowanego rzędu w zagadnieniach związanych z odlewnictwem 39 3.1. Zastosowanie metody KR w podzbiorach węzłów. 39. 3.2 Definicja kwadratury różniczkowej sterowanego rzędu (ang. Rank Controlled Differential Quadrature) 41 3.3 Budowa schematów rozwiązywania problemów transportu ciepła przy pomocy aproksymacji metodą KRSR 44 3.3.1. Dyskretyzacja pochodnej względem czasu. 3.3.2 Dyskretyzacja Fouriera - Kirchhoffa 3.3.3. pochodnych. przestrzennych. 44 w równaniu. różniczkowym 45. Dyskretyzacja warunków brzegowych przy pomocy metody KRSR. 47. 3.3.3.1. Warunek brzegowy I rodzaju. 47. 3.3.3.2. Warunek brzegowy II rodzaju. 48. 3.3.3.3. Warunek brzegowy III rodzaju. 49. 3.3.3.4. Warunek brzegowy IV rodzaju. 50 Strona 2.

(5) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. Rozdział 4: Analiza metodą KRSR. dokładności. przybliżenia. pochodnej 52. 4.1. Wpływ liczby węzłów siatki na dokładność przybliżenia. 53. 4.2. Analiza rozkładu błędu aproksymacji w dziedzinie. 58. 4.3. Analiza wpływu wartości rzędu minimalnego i maksymalnego metody KRSR. 60. Rozdział 5: Weryfikacja rozwiązania przybliżonego dla problemu transportu ciepła o znanym rozwiązaniu dokładnym 66 5.1. Model matematyczny przewodzenia ciepła w nieskończonej płycie. 66. 5.2. Rozwiązanie dokładne. 67. 5.3. Dyskretyzacja problemu. 68. 5.3.1. Metoda kwadratur różniczkowych sterowanego rzędu. 68. 5.3.2. Metoda różnic skończonych. 68. Wyniki modelowania numerycznego. 69. 5.4. Rozdział 6: Schematy krzepnięcia. umożliwiające. rozwiązywanie. 6.1 Model matematyczny krzepnięcia odlewu nieskończonej z wybranego stopu w formie piaskowej 6.2. modeli 75. płyty wykonanej 76. Numeryczny model krzepnięcia nieskończonej płyty. 78. 6.2.1. Dyskretyzacja pochodnych względem czasu oraz zmiennych przestrzennych 79. 6.2.2. Model krzepnięcia stopu, gdy znana jest funkcja fS. 6.2.2.1. 81. Algorytm obliczeniowy. 6.2.3 Dyskretna postać nieskończonej płyty. równań. 83 różniczkowych. opisujących. krzepnięcie 83. 6.3 Zastosowanie termicznej analizy różnicowej (DTA) do wyznaczenia przebiegu funkcji źródła ciepła przemian 84 6.3.1. Określenie szybkości przemiany. 85. 6.3.1.1 Teoretyczna analiza równań bilansu ciepła w celu wyprowadzenia zależności całkowej dla fS 85 6.3.1.2 Algorytm wyznaczania przebiegu funkcji udziału objętości zakrzepłej na podstawie danych w postaci dyskretnych wartości 86 6.4. Opis eksperymentu. 6.4.1. Metodyka badań eksperymentalnych. 87 87. 6.4.2 Wyznaczenie przebiegu funkcji źródła ciepła przemian na podstawie danych eksperymentalnych 88 Strona 3.

(6) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 6.4.3 Parametry symulacji - warunki jednoznaczności dla przypadku krzepnięcia nieskończonej płyty wykonanej ze stopu AZ91 w formie piaskowej 92 6.4.4. Wyniki modelowania numerycznego. Rozdział 7:. Dyskusja wyników badań. 93. 98. 7.1 Opracowania teoretyczne i eksperyment numeryczny podstawą optymalizacji parametrów metody KRSR 99 7.2. Zastosowanie metody KRSR podczas analizy procesów zależnych od czasu. 103. 7.3. Uwagi i wnioski końcowe. 109. Dodatek A: Aproksymacja parametrów termofizycznych dla stopu AZ91 oraz materiału formy 111 Funkcje przybliżające parametry termofizyczne stopu AZ91 wykorzystanego podczas obliczeń numerycznych 111 Ciepło właściwe. 111. Współczynnik przewodzenia ciepła. 112. Gęstość. 112. Funkcje przybliżające parametry termofizyczne dla materiału formy. 112. Ciepło właściwe. 112. Współczynnik przewodzenia ciepła. 112. Gęstość. 113. Literatura. 114. Streszczenie w języku polskim Streszczenie w języku angielskim. Strona 4.

(7) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. Rozdział 1 Podstawy teoretyczne kwadratury różniczkowej oraz matematycznego opisu procesów przewodzenia ciepła. Zrozumienie praw rządzących procesami fizycznymi jest coraz większe. Dziś prawie dla wszystkich problemów inżynierskich możliwe jest stworzenie modeli matematycznych, które bazując na prawach fizycznych opisują zmianę istotnych parametrów procesu. Model taki, opiera się na zestawieniu układu równań różniczkowych cząstkowych (RRCz) opisujących analizowane procesy, wraz z warunkami jednoznaczności (warunki geometryczne, fizyczne, brzegowe, dla problemów dynamicznych: początkowe). Przykładem takiego modelu może być równanie różniczkowe Fouriera-Kirchhoffa (FK), które uzupełnione o równania ciągłości masy, bilansu ciśnienia oraz stężeń z bardzo dobrym przybliżeniem opisują przepływ metalu we wnęce formy oraz zmianę jego temperatury. Ponadto jeżeli nieliniowy człon opisujący wewnętrzne źródło ciepła, zostanie wyrażony przez odpowiednio dobraną funkcję, opisującą kinetykę uwalniania energii przemiany, otrzymamy model, który pozwoli przewidywać mikrostrukturę wykonanego detalu. Niestety, znalezienie analitycznego rozwiązania tych równań jest często niemożliwe, a zawsze jest procesem bardzo trudnym. Zwykle metody znajdowania dokładnego rozwiązania można zastosować tylko dla kilku szczególnych przypadków. Jednak do opisu rzeczywistych procesów Strona 5.

(8) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. stosowane są układy RRCz o dużo większym stopniu skomplikowania. Ze względów praktycznych znalezienie rozwiązań problemów, z którymi borykają się inżynierowie jest bardzo istotne. Przewidywanie błędów technologii pozwala zakładom produkcyjnym zaoszczędzić znaczne kwoty na etapie planowania produkcji nowych elementów. Za zmniejszeniem ilości braków oraz oszczędnością materiałów zużytych do produkcji idzie zmniejszenie negatywnego wpływu zakładu na środowisko naturalne. Przewidywanie parametrów użytkowych na etapie symulacji przebiegu procesu pozwala określić czas bezpiecznego użytkowania detalu. Jedyną drogą pozostaje zatem poszukiwanie aproksymacji rozwiązania badanego modelu. Dlatego tak ważnym jest ciągłe rozwijanie metod przybliżonego rozwiązywania równań różniczkowych cząstkowych. Zwykle przybliżone rozwiązanie problemu, to zbiór ciągów wektorów o współrzędnych, które określają położenie analizowanych punktów w przestrzeni, związane z nimi wartości poszukiwanej funkcji oraz położenie na osi czasu (dla problemów dynamicznych). Punkty, w których znajdujemy przybliżone rozwiązanie nazywamy węzłami siatki. Na tym etapie dążenie do znalezienia rozwiązania modelu zapisanego przy pomocy układu równań różniczkowych cząstkowych wymaga postawienia i znalezienia odpowiedzi na pytanie: jaki związek występuje pomiędzy pochodną cząstkową, a poszukiwaną wartością funkcji w węźle siatki. Okazuje się, że istnieje takie powiązanie. Jest to numeryczna dyskretyzacja zagadnienia, pozwala ona na wykorzystanie informacji zawartych w pochodnych funkcji do znalezienia rozwiązania numerycznego [1]. Duży wpływ na dokładność oraz czas po jakim otrzymamy rozwiązanie postawionego zagadnienia ma wybór odpowiedniej metody numerycznej. Rozwój technologii informacyjnej pozwala na wykonanie symulacji wielu procesów zachodzących podczas zalewania oraz krzepnięcia odlewów. Rozważanie kilku procesów, które trwają równolegle i wpływają na siebie, stawia jednak wyższe niż dotychczas wymagania dotyczące dokładności prowadzonych obliczeń. Nawet niewielkie błędy pojawiające się w poszczególnych schematach kumulują się i zwielokrotniają w trakcie analizy procesów równoległych. Oznacza to, iż wciąż ważne jest poszukiwanie nowych metod numerycznych oraz analiza ich właściwości, a także określanie dla konkretnych metod obszarów zastosowań, w których są najskuteczniejsze. Prezentowana tutaj metoda kwadratur różniczkowych sterowanego rzędu (KRSR), rozszerzająca metodę kwadratur różniczkowych (KR), ze względu na swoje właściwości może stać się metodą odpowiednią do zastosowania w zagadnieniach związanych z transportem ciepła i masy. Metoda KRSR może być rozumiana jako uogólnienie metod różnic skończonych (RS), które to metody znalazły poczesne miejsce wśród przybliżonych metod stosowanych przez inżynierów. Duże znaczenie ma określenie parametrów, przy których zastosowanie metody KRSR da w efekcie istotną poprawę dokładności wyniku obliczeń numerycznych, nie wydłużając ich w nadmiernym stopniu. Znajomość tych parametrów pozwala indywidualnie potraktować każdy problem i konstruować dla niego odpowiedni schemat numeryczny. Dzięki temu uzyskane przybliżenie będzie najlepsze oraz zostanie uzyskane po dopuszczalnie długo trwających obliczeniach. Strona 6.

(9) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.1 Podstawy matematyczne opisu przewodzenia ciepła 1.1.1 Podstawowe równania rządzące procesami migracji ciepła Pod względem matematycznym przewodzenie ciepła w obszarze  opisane jest równaniem różniczkowym cząstkowym rzędu drugiego typu parabolicznego [2]. W zależności od postawionego problemu (obecność przemian lub zewnętrznych źródeł ciepła) równanie, to może być liniowe bądź nie. Wyprowadzenie równania różniczkowego przewodnictwa zwanego także równaniem energii lub równaniem Fouriera - Kirchhoffa można znaleźć w wielu pracach [3-6]. W tej pracy zostanie przytoczone równanie FK, w dwóch formach. Pierwsza postać bierze pod uwagę konwekcyjne ruchy medium powodowane różnicami gęstości podobszarów, które wywołane są gradientem temperatury:. cp . T  c p  u  div   grad T   qV , . (1.1). gdzie: cp [J kg-1 K-1] oznacza ciepło właściwe przy stałym ciśnieniu,  [W m-1 K-1] jest współczynnikiem przewodzenia ciepła,  [kg m-3] jest gęstością materiału, u [m s-1] jest wektorem prędkości ruchu medium w każdym kierunku układu, qV [W m-3] jest całkowitą wydajnością wewnętrznych źródeł ciepła. Postać ta (1.1) pozwala na opis pola temperatury w cieczach i gazach. Jeżeli mamy do czynienia z zagadnieniami, w których opisywana jest temperatura w ciele stałym, bądź ze względu na jej znikomy wpływ na przebieg zmian temperatury, można pominąć konwekcję, stosujemy uproszczoną postać równania FK:. cp . T  div   grad T   qV . . (1.2). Równanie Fouriera - Kiechhoffa w przedstawionej powyżej formie (1.2) opisuje pole temperatury w obszarze , w którym ciepło przenoszone jest tylko przez przewodzenie. Aby matematyczny opis przepływu ciepła w rozważanym obszarze  był kompletny, równanie FK musi zostać uzupełnione o tzw. warunki jednoznaczności. Wśród warunków jednoznaczności wyszczególnia się: warunki brzegowe, początkowe, geometryczne i fizyczne [3-6].. 1.1.2 Opis warunków jednoznaczności W paragrafie tym zostaną skrótowo przedstawione warunki niezbędne do pełnego opisu modelu, dla którego analizowany jest proces przewodzenia ciepła. Zamodelowanie takiego układu wymaga przede wszystkim doboru odpowiedniego prawa (lub praw) fizycznego opisującego badany proces. Następnie na bazie podstaw fizycznych wyprowadzane są równania, które opisują ilościowo zamianę ciepła w rozważanym. Strona 7.

(10) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. obszarze , równania: (1.1) bądź (1.2). Wspomniane równania muszą być wzbogacone o warunki jednoznaczności, by w pełni opisywać analizowany układ. 1.1.2.1 Warunki geometryczne Pojęcie to oznacza geometrię rozpatrywanego układu. W przypadku obszaru niejednorodnego podział na podobszary. Orientację obiektu oraz związanie go z konkretnym układem współrzędnych. 1.1.2.2 Warunki fizyczne Warunki fizyczne, to zbiór zależności określających wartości parametrów identyfikujących właściwości materiałów z jakich wykonane są kolejne podobszary. Zaliczamy do nich także parametry charakterystyczne dla zachowań danych parametrów. Wartości graniczne, punkty przejścia, itp. W przypadku przewodzenia problemów ciepła warunki fizyczne to przede wszystkim parametry termofizyczne podobszarów. Temperatury, przy których rozpoczynają bądź kończą się przemiany. Charakterystyczne dla rozważanego procesu temperatury. Jeżeli parametry fizyczne nie odpowiadają wartościom dla danego materiału, a są wynikiem modyfikacji związanych z przekształceniami wynikającymi z budowania modelu matematycznego ich opis. Taka sytuacja ma miejsce w przypadku stosowania metody zastępczej pojemności cieplnej podczas opisu krzepnięcia [3, 7]. 1.1.2.3 Warunki początkowe Warunki początkowe opisują pole poszukiwanego parametru we wszystkich podobszarach rozpatrywanego obszaru, w chwili, którą uznajemy za początkowy moment dla procesu. 1.1.2.4 Warunki brzegowe Równania te służą do opisu warunków jakie panują na brzegu obszaru oraz jakie panują na styku podobszarów stanowiących modelowany układ. Obszar  składa się z podobszarów i (. . i.  , i   j  i , j , i  j ), które charakteryzują się odmiennymi. i. właściwościami termofizycznymi lub otoczone są przez różne media. Dla obszarów tych rozdziela się myślowo brzegi obszaru : wyszczególnionych podobszarów ograniczających układ i oraz brzegi łączące niejednorodne podobszary i,j. Powstaje w ten sposób zbiór brzegów k, do których przypisuje się warunki brzegowe. W zależności od charakteru fizycznego brzegów w problemach przewodzenia ciepła wyróżnia się cztery rodzaje warunków brzegowych. W paragrafie (§3.3.3) zostaną przytoczone warunki brzegowe wyróżniane w odlewnictwie. Dla każdego typu zostanie przedstawiona metoda przybliżania ich metodą KRSR.. Strona 8.

(11) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.1.3 Równanie różniczkowe Fouriera - Kirchhoffa dla obszarów jednowymiarowych W niektórych przypadkach istnieje możliwość sprowadzenia rozważanego problemu przewodzenia ciepła do problemu jednowymiarowego. Taka sytuacja ma miejsce dla geometrii sprowadzalnych do przypadku nieskończonej płyty, nieskończonego walca oraz powłoki kulistej [3, 4]. Wzory te w postaci ogólnej mają postać: cp . T x,  1   m T x,    m  x   qV ,  x x  x . (1.3). gdzie: m = 0, 1, 2 odnoszą się odpowiednio do przypadków: płyty, walca i kuli [3], x [m] jest zmienną w przestrzeni ,  [s] oznacza czas. Zmienna x jest określona w pewnym przedziale [a, b]. Przedział ten jest dziedziną zadania i wielkość jego wynika z założeń modelu matematycznego przepływu ciepła. Dla przypadku płyty nieskończonej, m = 0, przy założeniu stałej wartości współczynnika przewodzenia ciepła równanie FK można zapisać w uproszczonej formie: T x,   2T x,  qV a  ,  x 2 cp . gdzie: a . (1.4).  [m2 s-1] jest współczynnikiem wyrównywania temperatury. cp . 1.1.3.1 Przypadek szczególny: równanie transportu ciepła bez funkcji źródła Często spotkać można się z procesami, w których nie zachodzi przemiana, zatem funkcja źródła występująca w równaniu FK jest równa zeru. W takim przypadku równanie różniczkowe przewodzenia ciepła znacząco się upraszcza. Element zawierający qV, stanowiący składnik równania decydujący o jego nieliniowości znika. Równanie przekształca się do postaci liniowej. W przypadku problemu przewodzenia ciepła w nieskończonej płycie równanie (1.3) przyjmuje postać: cp . T x,    T x,     ,  x  x . (1.5). T x,   2T x,  a .  x 2. (1.6). równanie (1.4) jest postaci:. W dalszej części niniejszej pracy zostaną przedstawione wyniki eksperymentów numerycznych. Obliczenia te prowadzone będą dla szczególnego przypadku równania FK (1.6), a także dla przypadku uwzględniającego krzepnięcie i zmienność parametrów Strona 9.

(12) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. termofizycznych (1.3). Posłużą one do prezentacji procedur umożliwiających wykorzystanie metody KRSR w zagadnieniach przewodzenia ciepła. Zestawienie wyników symulacji z rozwiązaniem dokładnym, dla specyficznego problemu dowiedzie przydatności analizowanej metody dla rozwiązywania zagadnień często spotykanych w odlewnictwie. W przypadku numerycznego modelowania krzepnięcia wyniki symulacji komputerowej z wykorzystaniem KRSR będą zestawione z wartościami zgromadzonymi podczas eksperymentów.. 1.1.4 Funkcja udziału objętości zakrzepłej Całkowitą wydajność wewnętrznych źródeł ciepła można związać ze zmianą udziału objętości zakrzepłej w badanym podobszarze o objętości V [m3]. Zakłada się, że w elemencie siatki o szerokości x [m], która odpowiada objętości V = F x po czasie  [s] od momentu rozpoczęcia obliczeń zakrzepła objętość metalu równa V S [m3]. Omawiane powyżej założenie prowadzi do równania [3]:. T x,     T x,    T x,  V S cp  V   F x  L ,  x  x  . (1.7). gdzie: L [J kg-1] jest utajonym ciepłem przemiany, F [m2] jest powierzchnią przez którą przechodzą linie strumienia ciepła. Po podzieleniu równania (1.7) przez V pojawia się w nim iloraz mówiący o tym jaka część całkowitej rozważanej objętości elementu V jest już zakrzepła: fS . V S . V. (1.8). Ze sposobu określenia funkcji fS wynika, że jest ona udziałem objętościowym objętości zakrzepłej w całej objętości rozpatrywanego podobszaru V, dlatego funkcję tę nazywa się funkcją udziału objętości zakrzepłej. Jasnym jest, że dla temperatury likwidus funkcja fS przyjmuje wartość 0, dla temperatury solidus fS osiąga 1. W zakresie krzepnięcia funkcja ta rośnie wraz ze spadkiem temperatury w rozpatrywanym elemencie. Po wykonaniu przejścia granicznego 0 i odpowiednim uproszczeniu V oraz F [3, 4, 7] równanie w formie jednowymiarowej dla przypadku nieskończonej płyty (1.7) przekształca się do postaci: cp . f x,  T x,    T x,   .    L S  x  x  . (1.9). Proces krzepnięcia odbywa się w przedziale temperatury, intuicyjnie jasnym jest, że funkcja fS musi być zależna od temperatury. Korzystając ze wzoru na pochodną funkcji złożonej otrzymuje się:. Strona 10.

(13) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. f S x,  f S T  T x,   .  T . (1.10). Równanie różniczkowe FK zapisane w postaci (1.9), w którym fS opisane jest związkiem (1.10) umożliwia stosowanie jednego z wielu algorytmów aproksymacji członu nieliniowego. Są to modele, które poprawiają lub bilansują temperaturę w rozważanym elemencie dyskretnym. Przykłady takich metod można odnaleźć w pracach [3, 7-10]. Funkcja fS (T) może zostać zdefiniowana a priori, uzyskana na podstawie danych empirycznych [4, 11] lub wyznaczona po zastosowaniu mikro modeli krzepnięcia [8-10, 12, 13].. 1.1.5 Zastępcza pojemność cieplna strefy dwufazowej Podejście prezentowane w tym podrozdziale prowadzi na drodze odpowiednich przekształceń do pojawienia się w równaniu energii nowego parametru termofizycznego. Zwany on jest zastępczą pojemnością cieplną strefy dwufazowej. Równocześnie równanie FK przyjmuje postać liniową. Występująca w równaniu (1.10) pochodna temperatury względem czasu może być powiązana z taką samą pochodną występującą po lewej stronie równości (1.9): f T   T x,    T x,       cp  L S  . T   x  x  . (1.11). Po takim przekształceniu równanie (1.9) przyjmuje formę "bezźródłową":. C. T x,    T x,     ,  x  x . (1.12). gdzie:. C  cp  L. f S T   J  , , T  kg K . (1.13). jest zastępczą pojemnością cieplną strefy dwufazowej. Podejście prezentowane powyżej prowadzi do metody przekształcania nieliniowego równania FK do postaci liniowej. Równanie FK zapisane w postaci (1.12) można rozwiązywać analogicznie jak równanie (1.5). Wpływ wewnętrznych źródeł ciepła na pole temperatury modyfikuje wartość parametru określonego jako zastępcza pojemność cieplna (1.13) co powoduje jego drastyczną zmianę w przedziale temperatur, w którym zachodzi krzepnięcie. Omówiona metoda linearyzacji równania różniczkowego Fouriera Kirchhoffa, mająca na celu umożliwienie zastosowania znacznie prostszych metod poszukiwania jego rozwiązania zwana jest metodą zastępczej pojemności cieplnej strefy dwufazowej.. Strona 11.

(14) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.1.6 Podsumowanie Przedstawione powyżej równania różniczkowe cząstkowe, w zależności od przebiegu modelowanego procesu przepływu ciepła tworzą układy równań. Ich rozwiązanie analityczne nawet w bardzo prostych przypadkach nie jest zwykle możliwe. W takich zagadnieniach stosuje się schematy bazujące na numerycznych metodach aproksymacji rozwiązania zadanego układu RRCz, te metody to między innymi: metoda różnic skończonych (MRS) [7, 14-16], metoda elementów skończonych (MES) [17, 18], metoda elementów brzegowych (MEB) [19, 20], metoda objętości skończonych (MOS) [21]. W niniejszej rozprawie zostanie przedstawiona metoda kwadratur różniczkowych (KR) oraz jej autorska modyfikacja zwana metodą kwadratur różniczkowych sterowanego rzędu (KRSR), a także sposoby w jaki można ją stosować podczas rozwiązywania różnych zagadnień opisujących transport ciepła.. 1.2 Opis metody kwadratur różniczkowych Metoda kwadratur różniczkowych została po raz pierwszy zaprezentowana przez Richarda Bellmana i wsp. [22-24] na początku lat '70 XX wieku. Jest to metoda dyskretyzacji pochodnych w równaniach różniczkowych. Jej idea jest zbliżona do idei metody całkowania numerycznego zwanej metodą kwadratur całkowych [17, 18, 25], w których całkę oznaczoną pewnej funkcji, przybliża się kombinacją liniową wartości tej funkcji, w ustalonych węzłach siatki. Metoda KR opiera się na podobnym pomyśle: pochodną funkcji w punkcie siatki zastępuje się sumą ważoną wartości poszukiwanej funkcji, we wcześniej wybranych punktach dyskretyzacji. Bellman wraz ze współpracownikami [22, 23] zaproponował dwie metody poszukiwania współczynników wagowych metody KR. Pierwsza polegała na rozwiązywaniu układu równań algebraicznych. Zaletą tego podejścia jest to, że można je stosować dla dowolnie wybranych węzłów siatki. Wadą natomiast to, że znalezienie ich wymaga rozwiązania układu równań z macierzą Vandermonde’a, która jest bardzo źle uwarunkowana [25]. Druga z podanych przez Bellmana metod przedstawia bardzo prostą formułę znajdowania współczynników wagowych opartą na wielomianach Legendre’a [1, 23]. Wadą tej metody jest jednak to, że węzły siatki są zdeterminowane przez pierwiastki wielomianu Legendre’a określonego w odpowiednim przedziale. W związku z trudnościami przy wyznaczaniu współczynników wagowych metoda KR była stosowana sporadycznie. Niemal w trzydzieści lat po opracowaniu podstaw metody KR przez Bellmana, Quan i Chang [26, 27] zaproponowali kolejną metodę wyznaczania współczynników wagowych metody, była ona oparta na wielomianach interpolacyjnych Lagrange’a. Pozwoliła ona na opracowanie jawnych wzorów opisujących wartości współczynników metody KR dla pochodnych pierwszego i drugiego rzędu. Kolejny przełom w rozwoju metody KR stanowiły prace Shu i Richards’a [28-31], w których autorzy uogólniają wyznaczanie współczynników wagowych dzięki zastosowaniu teorii aproksymacji funkcji w przestrzeniach wielomianów wysokiego stopnia. Podejście to pozwoliło na opracowanie jawnych wzorów rekurencyjnych pozwalających wyznaczyć współczynniki kwadratury dla pochodnej dowolnego rzędu. Inaczej niż w podejściu Bellmana wzory te są prawdziwe dla dowolnych rozkładów węzłów siatki. Strona 12.

(15) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.2.1 Definicja kwadratury różniczkowej Niech f :   a, b  x  f x    będzie funkcją klasy. C. m+1. . Na odcinku [a, b]. wprowadzamy N węzłową siatkę: S   xi : a  x1  x2  x3    xN 1  xN  b .. (1.14). DEFINICJA 1 Kwadraturą różniczkową (KR) nazywa się metodę przybliżania pochodnych funkcji f zdefiniowaną następującymi wzorami:. xi    i(,mj) f x j  , N. f. ( m). i, j = 1, 2, 3,..., N,. (1.15). j 1. gdzie:. f ( m) xi  oznacza m-tą pochodną funkcji f względem zmiennej x, w punkcie xi,. i(,mj ) są współczynnikami wagowymi metody KR dla m-tej pochodnej. Wzór (1.15) można zapisać w postaci macierzowej:. 1(,m1 ) 1(,m2 )  1(,mN)   f x1    f ( m ) x1    (m)    ( m) 2( m, N)   f x2    f ( m ) x2   2,1 1, 2             .    i(,mj )                   (m)      (m)   N( m, N)   f xN   f ( m ) xN   N ,1  N , 2. (1.16). Zapisana w postaci macierzowej definicja kwadratury różniczkowej w sposób bardzo intuicyjny obrazuje, na czym polega zastosowanie metody KR podczas rozwiązywania RRCz. Pochodną w punkcie xi występującą w równaniu różniczkowym zastępuje się iloczynem skalarnym wektora utworzonego z i-tego wiersza macierzy współczynników wagowych przez wektor wartości funkcji we wszystkich węzłach siatki S. W dalszej części pracy macierz współczynników wagowych dla metody KR zastosowanej do pochodnej rzędu m oznaczona będzie symbolem [W(m)]. Kluczowym zagadnieniem związanym z przybliżonym rozwiązywaniem RRCz z zastosowaniem metody KR jest znalezienie współczynników wagowych i(,mj ) kwadratury.. Strona 13.

(16) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.2.2 Wyznaczanie współczynników wagowych metody KR dla pochodnej pierwszego rzędu Bardzo istotnym zagadnieniem podczas stosowania metody KR jest opracowanie skutecznego sposobu wyznaczania jej współczynników wagowych. Definicja kwadratury różniczkowej (Definicja 1) może stanowić podstawę algorytmu poszukiwania ich wartości. Wybrany układ wielomianów, będący bazą przestrzeni VN, może stanowić układ funkcji testowych, dla których przeprowadzana jest aproksymacja metodą KR. Współczynniki wagowe kwadratury chcemy dobrać tak, by była ona dokładna dla każdej funkcji bazowej i dla każdego węzła siatki S. Przeprowadzenie procedury opisanej definicją (1.15) dla wszystkich funkcji bazowych prowadzi do układu N × N niezależnych równań algebraicznych. Rozwiązanie tego układu równań pozwala określić wartości współczynników wagowych kwadratury różniczkowej. 1.2.2.1 Baza kanoniczna przestrzeni wielomianów Bellman i wsp. [22] zaproponowali podejście polegające na wykorzystaniu jako układu funkcji testowych wielomianów bazy kanonicznej przestrzeni wielomianów stopnia naturalnego mniejszego niż N. Niech funkcje bazy kanonicznej będą oznaczone przez:. ek  x k 1 , k  1, 2, 3,, N .. (1.17). Podstawiając te funkcje do równania (1.15) otrzymujemy następujący układ równań algebraicznych:.  N (1)  i , j  0,  j 1  N (1) dla i  1, 2, 3,  , N .  i , j x j  1,  j 1 N  i(,1j) x kj  kxi( k 1) , k  2, 3,  , N  1,  j 1. (1.18). Przedstawiony powyżej układ równań można rozwiązywać po kolei dla kolejnych i, otrzymując w ten sposób kolejne wiersze macierzy współczynników metody KR. [W(1)] (1.16).. Jak to było już wspomniane, układy równań (1.18) są układami równań z macierzą Vandermonde’a, która jest bardzo źle uwarunkowana. W przypadku tym nawet zaburzenia dokładności na poziomie reprezentacji danych mogą prowadzić do uzyskania. . błędnych wartości współczynników i(,1j). . i , j 1, 2,...,N. . Wpływ złego uwarunkowania macierzy. rośnie bardzo znacznie wraz ze wzrostem liczby węzłów siatki dyskretnej N. Podejście to nie sprawdza się zatem w przypadku zastosowań praktycznych. Jednak, zastosowanie go w przypadku problemów nie wymagających gęstych siatek pokazało bardzo obiecujące własności metody KR [23, 32, 33]. Zwróciło to uwagę osób zajmujących się Strona 14.

(17) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. rozwijaniem metod numerycznych na metodę kwadratur różniczkowych co skutkowało propozycjami jej ulepszenia (§1.2.2.3). Opierały się one na analizach efektów stosowania układów funkcji różnych typów jako zbiorów funkcji testowych. Możliwość uzyskania współczynników wagowych metody kwadratur różniczkowych przy użyciu innej bazy wielomianów stopnia mniejszego niż N prowadzi do pytania o jednoznaczność wartości współczynników wagowych. 1.2.2.2 Jednoznaczność współczynników wagowych KR dla różnych baz [1] Niech VN będzie przestrzenią wektorową wielomianów stopnia mniejszego lub równego N, zbiory B1  1 x ,  2 x ,,  N x  i B2  1 x ,  2 x ,,  N x  będą różnymi bazami wielomianów w tej przestrzeni. Niech [Q] oznacza macierz przejścia z bazy, B1 do B2. Zastosowanie definicji KR (1.16) dla wielomianów baz B1, B2 prowadzi do wzorów:.  B x  B S W ,  B x  B S W , ( m) 1. k. 1. (k ) 1. dla k  1, 2,..., N ,. (1.19). ( m) 2. k. 2. (k ) 2. dla k  1, 2,..., N ,. (1.20). gdzie:. W  – oznacza. -tą kolumnę macierzy współczynników wagowych dla pochodnej. (k ). i. n-tego rzędu [W(n)] wyznaczoną przy użyciu wielomianów bazy Bi;. B1 S    i x j i, j1,2,...,N , B2 S   i x j i, j1,2,...,N. - są macierzami, w których w i-tym. wierszu i w j-ej kolumnie znajduje się wartość i-tego wielomianu bazowego w j-ym węźle siatki S.. B x   x ,  x ,,  x  ( m) 1. k. ( m) 1. ( m) 2. k. ( m) N. k. T. B x   x ,  x ,,  x  ( m) 2. k. ( m) 1. ( m) 2. k. k. ( m) N. ,. k. T. k. – są wektorami wartości pochodnych. m - tego rzędu funkcji bazowych w wybranym punkcie xk. Mnożąc równania (1.19) lewostronnie przez macierz przejścia [Q] otrzymujemy:. Q B1(m) xk  QB1 S W1(k ) ,. dla k  1, 2,..., N ,. jednak macierz przejścia ma następujące własności:. QB1(m)   B2(m) , stąd:.  B x  B S W , ( m) 2. k. 2. (k ) 1. (1.21). Q B1(m) xk   B2(m) xk .  k  1, 2,..., N .. oraz. (1.22). Wyrażenie (1.22) zestawione z (1.19) daje wniosek: Strona 15.

(18) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. W1   W2 .. (1.23). Tym samym została udowodniona prawdziwość następującego twierdzenia: TWIERDZENIE 1 Współczynniki wagowe kwadratury różniczkowej dla aproksymacji m-tej pochodnej, i(,mj ) , nie zależą od wyboru bazy w przestrzeni wielomianów VN.. Twierdzenie to pozwala na sformułowanie wniosku, że wartości współczynników wagowych kwadratury będą takie same, bez znaczenia w jakiej bazie były one wyznaczone. Wniosek ten stał się kanwą poszukiwania nowych algorytmów wyznaczania współczynników wagowych metody KR. Opierały się one przede wszystkim na analizie skuteczności zastosowań różnych baz w przestrzeni wielomianów wysokiego rzędu do wyprowadzania wzorów opisujących wartości tych współczynników. 1.2.2.3 Baza wielomianów interpolacyjnych Lagrange’a Zdecydowanie najlepszy algorytm znajdowania współczynników wagowych KR, w kategoriach zastosowań praktycznych, przestawiony został w 1989 roku przez Quana i Changa [26, 27]. Polega on na wykorzystaniu jako funkcji testowych wektorów bazy wielomianów interpolacyjnych Lagrange’a: N. lk ( x)   i 1 ik. x  xi , xk  xi. (1.24). gdzie: xi należą do N węzłowej siatki S zdefiniowanej jak w (1.14), k = 1, 2, 3, ..., N. Wielomiany interpolacyjne mają następującą własność:. 1, gdy k  s, lk xs    0, gdy k  s.. (1.25). Wykorzystanie tych wielomianów jako funkcji testowych, sprowadza układ równań (1.15) do postaci:. i(,1k)  lk(1) xi , dla i, k  1, 2, 3,..., N .. (1.26). Pierwsza pochodna funkcji interpolacyjnej Lagrange’a wynosi:. Strona 16.

(19) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów.  N        x  xn  j 1  n 1  j k  nk , j  (1) lk ( x)  . N  xk  xi  N. (1.27). i 1 ik. Wyznaczanie wartości pierwszej pochodnej wielomianów interpolacyjnych w konkretnych punktach siatki, sprowadza to wyrażenie, do jednej z dwóch postaci, zależnej od indeksu węzła oraz indeksu wielomianu interpolacyjnego Lagrange’a:  przypadek, gdy indeks węzła jest taki jak indeks funkcji bazowej: N. k(1,k)   j 1 j k. 1 , xk  x j . (1.28).  przypadek, gdy indeksy węzła i funkcji bazowej różnią się:. i(,1k) . N x  x  1 i j .  xk  xi  j1 xk  x j . (1.29). j i , k. Powyższe wzory (1.28) i (1.29) tworzą algorytm, który pozwala wyznaczyć współczynniki kwadratury bez konieczności rozwiązywania skomplikowanego układu równań. Procedurę wyznaczania współczynników wagowych można usprawnić korzystając ze wzoru (1.18) dla e0 (x) oraz twierdzenia 1: N  (1) xi  xn  , dla i  j 1    i , j  x j  xi  n1 xk  xn   n i ,k  N  (1)    (1) .  i, j  i ,i j 1  j i. (1.30). Na podkreślenie zasługuje fakt, iż wartości współczynników wagowych metody KR zależą tylko od współrzędnych siatki S wprowadzonej w dziedzinie zadania. Macierz współczynników metody KR jest centro-symetryczna. Można wykorzystać tę własność do redukcji operacji arytmetycznych koniecznych do jej wyznaczenia. Algorytm taki został opisany przez Chen’a i wsp. [34]. Jest on bardzo przydatny w przypadku problemów, w których podczas rozwiązywania zmienia się rozkład węzłów dyskretyzacji. Co może być związane z dynamiką samego procesu (np. problemy o ruchomej siatce) własność centro-symetryczności skraca wówczas czas konieczny do wyznaczenia współczynników wagowych.. Strona 17.

(20) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.2.3 Wyznaczanie współczynników metody KR dla pochodnych drugiego i wyższych rzędów Wielomiany interpolacyjne Lagrange’a znajdują również zastosowanie w przypadku dyskretyzacji pochodnych drugiego i wyższych rzędów. Podejście takie pozwala otrzymać wzory rekurencyjne na współczynniki dyskretyzacji pochodnych kolejnych rzędów. Możliwe jest również zastosowanie własności pochodnej i wyznaczanie macierzy współczynników wagowych wyższych rzędów poprzez iteracyjne wykonywanie odpowiedniej ilości iloczynów macierzy współczynników (1.16) dla pochodnej rzędu pierwszego. 1.2.3.1 Metoda oparta o własności macierzy współczynników wagowych Operacja wyznaczania kolejnych pochodnych cząstkowych przy pomocy pochodnej cząstkowej pierwszego rzędu, może być zapisana rekurencyjnie:. m f    m1 f   x m x  x m1.  , dla m  2, 3,... . . (1.31). Metoda aproksymacji pochodnej metodą KR pozwala pochodne po obu stronach równości zastąpić kombinacjami liniowymi wartości funkcji f lub jej odpowiedniej pochodnej, w węzłach siatki i odpowiednich współczynników wagowych: f ( m ) xi    i(,mj ) f x j   N. j 1 N.   k 1. N. (1) i ,k.  j 1. N  ( m1) xi    i(,1k) f ( m1) xk   f x k 1. ( m1) k, j. f x j    N. j 1. N.  k 1. . (1) i ,k. ( m 1) k, j. f x j ,. (1.32). stąd możemy wnioskować, że:. i(,mj)   i(,1k)k( m, j1)  W ( m)   W (1)  W m1   W (1)  m N. (1.33). k 1. Wyznaczanie współczynników wagowych przy użyciu wzoru (1.33) jest bardzo proste w implementacji numerycznej, jednak wymaga dużego nakładu obliczeń. Dlatego metoda ta bywa stosowana w przypadku siatek o niewielkiej liczbie węzłów oraz w niektórych algorytmach wprowadzania warunków brzegowych [35-37]. 1.2.3.2 Metoda rekurencyjnego wyznaczania współczynników wagowych dla pochodnych drugiego i wyższych rzędów Wielomiany interpolacyjne Lagrange’a zastosowane jako funkcje testowe w definicji KR (1.15) pozwalają wyznaczyć wzory rekurencyjne do wyznaczania współczynników wagowych aproksymacji pochodnych rzędów wyższych lub równych dwa. Podstawienie tych wielomianów do wzoru (1.15) i wykorzystanie własności (1.25) prowadzi do zależności:. i(,mk)  lk( m) xi , dla i, k  1, 2, 3,..., N ,. (1.34). gdzie: m oznacza rząd przybliżanej pochodnej. Strona 18.

(21) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. Podobnie jak w przypadku pochodnej pierwszego stopnia (§1.2.2.3) można wziąć pod uwagę dwa przypadki: przypadek, kiedy indeks węzła jest taki sam jak indeks wielomianu bazowego oraz przypadek, kiedy indeksy wielomianu interpolacyjnego i węzła są różne. Wzory rekurencyjne na wyznaczanie współczynników wagowych zostały przedstawione po raz pierwszy przez Quana i Changa [26, 27]. W przypadku, gdy indeks funkcji bazowej jest taki sam jak indeks węzła (diagonala macierzy współczynników), formuła określająca wartość współczynnika jest postaci:. i(,mi) . M m1 xi  , m  1 M (1) xi . (1.35). gdzie: M x   x  x1 x  x2    x  xN  .. (1.36). Ze względu na trudności z wyznaczeniem M (m+1)(x) wzoru tego nie stosuje się w zastosowaniach praktycznych. Podczas wyliczania wartości współczynników wagowych znajdujących się na przekątnej głównej korzysta się z własności, którą otrzymuje się po wstawieniu do definicji KR (1.15) funkcji em-1 (x), zdefiniowanej w (1.17), jako funkcji testowej: N.  j 1. ( m) i, j.  0.. (1.37). Współczynniki wagowe metody KR, i(,mj ) dla i ≠ j można związać rekurencyjnie ze współczynnikami i(,1j). oraz i(,mj1). [1, 26, 27]. Algorytm pozwalający rekurencyjnie. wyznaczyć współczynniki wagowe metody KR dla pochodnej rzędu m, uwzględniający właściwość (1.37) jest postaci:  ( m)  1 m1 i,mj1   i , j i ,i  , dla i  j   m  i, j     x  x  i j    N  ( m )    ( m ) . i, j  i ,i j 1  j i. (1.38). Algorytm opisany wzorem (1.38) wymaga znajomości wartości współczynników i(,1j) , które można otrzymać przy użyciu formuł (1.30). Wyznaczanie współczynników wagowych metody KR zgodnie ze wzorem (1.38) wymaga znacznie mniejszego nakładu obliczeń, niż w przypadku algorytmu (1.33). Współczynniki wagowe pochodnej rzędu m tak jak dla pochodnej pierwszego rzędu zależą wyłącznie od współrzędnych węzłów siatki S.. Strona 19.

(22) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. 1.3 Dokładność przybliżenia metodą KR Pełna teoretyczna analiza czynników wpływających na dokładność metod numerycznych opartych kwadraturach różniczkowych nie została dotychczas przeprowadzona. Głównym źródłem informacji o charakterze pojawiających się błędów i wskazówek, w jaki sposób ich uniknąć są wyniki eksperymentów numerycznych przeprowadzanych bezpośrednio dla metody KR. W niniejszej rozprawie przedstawiono wybór istotnych twierdzeń i oszacowań, które stanowić mogą przyczynek do analizy teoretycznej błędu pojawiającego się podczas stosowania metody KR. Na podstawie wykazanych związków KR z różniczkowaniem wielomianu interpolacyjnego autor tej pracy dowodzi, iż nie jest możliwe uzyskanie aproksymacji wysokiej dokładności przy zastosowaniu metody KR dla licznych siatek dyskretnych. Wynik ten staje się jednak podstawą do sformułowania metody kwadratur różniczkowych sterowanego rzędu, która charakteryzuje się wysoką dokładnością nawet dla bardzo gęstych siatek.. 1.3.1 Zagadnienie najlepszej aproksymacji W paragrafie tym, zostaną przytoczone pewne, istotne ze względu na dalszy opis metody KR, pojęcia i twierdzenia dotyczące teorii aproksymacji. W dalszej części rozpatrywane będą zagadnienia związane z aproksymacją w unormowanej przestrzeni wektorowej funkcji rzeczywistych. Przedstawione tematy są bardzo istotne w kontekście analizy możliwości wpływania na dokładność numerycznego rozwiązywania problemów matematyczno – fizycznych. 1.3.1.1 Zagadnienie najlepszej aproksymacji w przestrzeni wektorowej Zadanie aproksymacji pojawia się, gdy konieczne jest wyrażenie pewnej funkcji prostszą (np. wyliczając wartości funkcji elementarnych przy pomocy kalkulatora, liczymy ich przybliżone sumy częściowe rozwinięć w szereg Maclaurina) lub poszukiwanie przybliżeń pewnych funkcji na podstawie informacji o jej wartościach lub wartościach pochodnych (wszelkie interpolacje oraz numeryczne rozwiązywanie RRCz). Zadanie aproksymacji liniowej zdefiniowane jest następująco [38]: DEFINICJA 2 Dla ustalonego elementu f z przestrzeni wektorowej unormowanej V szukamy elementu h*, należącego do danej podprzestrzeni wektorowej VN, skończenie wymiarowej, takiego że:. f  h*  f  h. (1.39). dla wszystkich h VN.. Błędem aproksymacji elementu f względem podprzestrzeni VN nazywamy wielkość VN ( f ) równą Strona 20.

(23) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów.  V  f   inf hU f  h ,. (1.40). N. a element h*, dla którego zachodzi równość f  h *   VN  f  elementem optymalnym. Poniższe twierdzenie pozwala na wprowadzenie pojęcia najlepszej aproksymacji: TWIERDZENIE 2 Niechaj V będzie unormowaną przestrzenią wektorową z normą.  , elementy. 1 , 2 ,,  N stanowią układ niezależnych liniowo elementów przestrzeni V. Elementy te rozpinają podprzestrzeń liniową VN  V. Dla pewnego f  V istnieje wówczas element N. VN  h*   ai*i , że: i 1. h VN : f  h *  f  h .. (1.41). Element h* jest najbliższym do elementu f spośród wektorów podprzestrzeni VN w sensie normy  . Element ten nazywamy elementem najlepszej aproksymacji. W zgodzie z założeniami twierdzenia 2 układ elementów i stanowi bazę przestrzeni wektorowej VN. Pojawia się teraz pytanie o jednoznaczność wyboru elementu najlepiej przybliżającego f. Odpowiedź na nie daje następujące twierdzenie: TWIERDZENIE 3 Niech VN będzie przestrzenią unitarną VN ,  generowana. iloczynem. skalarnym:. x . i norma w przestrzeni VN niech będzie. x, x .. Wówczas. element. najlepszej. aproksymacji, dla elementu f  VN, jest jedyny i określony następującą tożsamością:. h VN : f  h*, h  0 .. (1.42). Z punktu widzenia zastosowań w analizie numerycznej, nasuwa się wniosek, że zadanie aproksymacji nad pewną podprzestrzenią liniową prowadzi do wyznaczenia elementu najlepszej aproksymacji. 1.3.1.2 Typy aproksymacji Powyższe twierdzenie (Tw. 3) pozwala na konkluzję, iż wybór metody aproksymacji jest równoważny wyborowi metody pomiaru odległości w danej przestrzeni unitarnej. Różne podejścia do tego zagadnienia pozwalają na sformułowanie trzech, najczęściej wykorzystywanych w praktyce, metod aproksymacji:. Strona 21.

(24) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów.  Aproksymacja jednostajna; polega na poszukiwaniu funkcji, której maksymalne odchylenie od aproksymowanej funkcji będzie najmniejsze.  Aproksymacja średniokwadratowa; polega na poszukiwaniu funkcji, której suma kwadratów odchyleń od aproksymowanej funkcji, w wybranych punktach będzie najmniejsza.  Interpolacja; polega na poszukiwaniu funkcji, która dokładnie przybliża aproksymowaną funkcję w wybranych punktach. Można oczywiście określić inne metody aproksymacji. W dalszej części nacisk zostanie położony na interpolację wielomianową, gdyż gra ona kluczową rolę jako podstawa metody KR.. 1.3.2 Interpolacja wielomianowa Rozważamy przedział a, b   i pewną funkcję rzeczywistą f : a, b   , w przedziale tym, znajduje się N argumentów, xi, rozmieszczonych w dziedzinie zgodnie z (1.14), nazywanych dalej węzłami interpolacji. W węzłach tych znamy wartości funkcji i są one równe odpowiednio: f xi   f i , i  1, 2,, N . Zadanie interpolacyjne polega na znalezieniu takich parametrów a1 , a2 , , aN   , takich by we wszystkich węzłach interpolacji spełniona była równość: N. k  1, 2, , N :.  a  (x )  f i 1. i. i. k. k. ,. (1.43). gdzie: wektory i stanowią bazę przestrzeni wektorowej VN, w której poszukujemy elementu najlepszej aproksymacji. Typ interpolacji związany jest z wyborem bazy przestrzeni wektorowej VN. Podstawą dla obecnie stosowanych metod KR może być interpolacja trygonometryczna [39, 40], interpolacja wielomianami oraz interpolacja funkcjami sklejanymi [41, 42]. W tej pracy zostanie przedstawiona wyłącznie interpolacja wielomianami, a w dalszej części, metoda KR oparta na bazie wielomianów interpolacyjnych. Przyjmując jako funkcje aproksymujące wielomiany stopnia mniejszego niż N otrzymujemy zadanie interpolacji wielomianami. Interpolacja wielomianowa jest operacją, która rzutuje funkcję f (x), ciągłą na przedziale [a, b] na podprzestrzeń wielomianów rzędu niższego niż N, VN. Operacja interpolacji funkcji f, oparta na węzłach siatki S (1.14) oznaczone zostanie symbolem I NS  f  , jest to jest liniowa projekcja na podprzestrzeń VN. Jej bazą mogą być różnie dobrane wielomiany. Istnieją trzy podstawowe algorytmy znajdowania współczynników ai we wzorze (1.43), są to: algorytm klasyczny, algorytm Newtona oraz algorytm Lagrange’a [17, 18, 25]. Algorytmy interpolacji w przestrzeniach wielomianów skończonego stopnia opierają się na określonych układach wielomianów stanowiących bazy danej przestrzeni. Znalezione współczynniki interpolacji mogą różnić się co do wartości, gdyż Strona 22.

(25) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. są to współrzędne wektora w różnych bazach, jednak przy spełnieniu założeń poniższego twierdzenia wielomian interpolacyjny jest określony jednoznacznie. TWIERDZENIE 4 Istnieje dokładnie jeden wielomian stopnia N-1 przyjmujący określone wartości w N różnych węzłach interpolacji.. Twierdzenie to pozwala na wybór dowolnej metody interpolacji, gdyż otrzymany element najlepszej aproksymacji w przypadku N par danych jest jedyny. Wzór interpolacyjny Lagrange’a jest postaci: N. p N ( x)  I NS ( f )( x)   f ( xk )lk ( x) ,. (1.44). k 1. gdzie: lk (x) jest k-tym wielomianem interpolacyjnym Lagrange’a wyrażonym wzorem (1.24). Wielomiany interpolacyjne Lagrange’a zdefiniowane wzorem (1.24) stanowią bazę przestrzeni VN. 1.3.2.1 Dokładność interpolacji Twierdzenie Weierstrassa sugeruje możliwość wskazania takiego ciągu wielomianów, który dokładnie przybliża dowolną funkcję rzeczywistą f. TWIERDZENIE 5(Twierdzenie Weierstrassa o aproksymacji wielomianami) Niechaj f (x) będzie ciągłą funkcją określoną na domkniętym przedziale [a, b] o wartościach rzeczywistych. Można wówczas wskazać ciąg wielomianów stopnia N-1: pN (x), który jest jednostajnie zbieżny do funkcji f (x) na przedziale [a, b].. Twierdzenie Weierstrassa może również zostać wypowiedziane w inny sposób: jeżeli funkcja f (x) jest ciągłą funkcją o wartościach rzeczywistych, określoną na przedziale [a, b], to dla każdego  > 0, istnieje wielomian pN (x), o pewnym, określonym stopniu (N-1), N=N (), dla którego spełniona jest nierówność:. max f ( x)  pN ( x)   .. x[ a ,b ]. (1.45). W praktyce twierdzenie to, postuluje istnienie wielomianu wysokiego rzędu, który z dowolną dokładnością przybliża odpowiednio gładką funkcję f. Można zatem zapisać:. Strona 23.

(26) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów N. f x   p N x    ci x i 1 ,. (1.46). i 1. Wstawiając kolejne wartości xk i fk otrzymujemy układ równań algebraicznych z macierzą Vandermonde’a, który definiuje klasyczny algorytm interpolacji: N. k  1, 2, , N : p N xk    ci xki 1 .. (1.47). i 1. Układ równań z macierzą Vandermonde’a posiada dokładnie jedno rozwiązanie. Jak to już zostało wspomniane, jego rozwiązanie metodami numerycznymi przysparza jednak kłopotów, ze względu na złe uwarunkowanie zadania. Biorąc pod uwagę twierdzenie o jednoznaczności funkcji interpolującej (Tw. 4) zadanie to można rozwiązać korzystając ze wzorów explicite (1.44), (1.24). Nie znaczy to jednak, że można wykorzystać wielomiany interpolacyjne do znajdowania dowolnej dokładności przybliżenia zadanej funkcji f (x) dla x należących do przedziału [a, b]. Niech f. oznacza normę przestrzeni funkcji, zadaną wzorem:. f  max f ( x). (1.48). x[ a ,b ]. W przypadku, gdy funkcja f jest N razy różniczkowalna, a jej N-ta pochodna jest ciągła i określona w dziedzinie, prawdziwe jest oszacowanie (§1.3.3.1):. Mf. f ( x)  p N ( x) . N!. N.  (x  x ) , i. (1.49). i 1. gdzie: M f  f ( N )  max f ( N ) ( x) . x[ a ,b ]. Wartość wielomianu interpolacyjnego jest zależna nie tylko od wartości N-tej pochodnej interpolowanej funkcji w rozpatrywanym przedziale, ale także od liczby węzłów interpolacji, N. N, oraz zachowania wielomianu M x    x  xi  w przedziale [a, b]. Dlatego analiza i 1. zbieżności wielomianu interpolacyjnego jest złożonym zagadnieniem. Znaczącą rolę podczas analizy zbieżności wielomianów interpolacyjnych odgrywa funkcja Lebesgue’a [43-45]: N. SN ( x)  max I NS  f ( x)   li ( x) f 1. (1.50). i 1. oraz stała Lebesgue’a [46] (zwana też liczbą Lebesgue’a dla interpolacji), która ma wartość równą maksimum funkcji Lebesgue’a dla ustalonej siatki S:. SN I   SN  max SN x  . x[ a ,b ]. (1.51) Strona 24.

(27) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. Niech f * = f - p*, wówczas prawdziwe jest następujące oszacowanie:. max I NS  f *( x)  max xa ,b . xa ,b. N. f l x   k 1. * k k. N   f *  max  lk x   xa ,b  k 1.    .    f *  max  lk x    f * max SN x   SN I   E N  f  , xa ,b   xa ,b  k 1  N. (1.52). max I NS  f ( x)  f ( x)  max I NS  f ( x)  p * ( x)  p * ( x)  f ( x)  xa ,b .  max I xa ,b . gdzie:. xa ,b . S N. f.  . . . (1.53). .  p *( x)  E N  f   max SN ( x)  1 E N  f   SN ( I )  1 E N  f , xa ,b. EN  f   inf  max f  p   f  p * , natomiast pVN  xa ,b  . p * jest elementem najlepszej. aproksymacji w przestrzeni wielomianów stopnia mniejszego niż N. Powyższe oszacowania pozwalają na stwierdzenie, iż ciąg wielomianów interpolacyjnych. . . S  0 . I NS  f  zmierza jednostajnie do f (x), wtedy i tylko wtedy, gdy  N I   1 EN  f   N . Zgodnie z Twierdzeniem Banacha – Steinhausa jest to możliwe tylko w przypadku, gdy.  I   1 są ograniczone. Warunek ten nie może być jednak spełniony, skoro (jak pokazał S N. Faber G. [47]) dla dowolnego rozkładu węzłów siatki S w przedziale [-1, 1], istnieje dodatnia stała c:. SN ( I )  2 ln( N )  c .. (1.54). Metody dokładniejszej analizy wartości stałej Lebesgue’a były [45 ,48] i wciąż są poszukiwane. Obecnie za najdokładniejsze przyjmowane jest oszacowanie podane przez Vértesi'ego dla interpolacji w przedziale [-1, 1] [49]:. . . min SN ( I )  2 ln( N )  2 (  ln 4 )  O ( lnlnlnNN ) 2 , S. (1.55). gdzie:  jest stałą Eulera - Mascheroniego ( = 0,577215665..) [50], zatem człon jest w przybliżeniu równy 0,521251...,. O jest symbolem Landaua. 2 . (  ln 4 ). (mówi. się, że. df. f n  Og n   r , c   : n  c  f n  r g (n) ). Analizy przedstawione w tym paragrafie przeprowadzone były dla siatek złożonych z węzłów rozmieszczonych w przedziale [-1, 1]. Dla siatki o węzłach znajdujących się w dowolnym podzbiorze liczb rzeczywistych istnieje odwzorowanie bijektywne na siatkę o węzłach z przedziału [-1, 1]. Wielkość przedziału nie wpływa zatem na charakter zbieżności, zmiana długości przedziału zmodyfikuje natomiast wartości stałych występujących w prezentowanych formułach. Przedstawione powyżej rozumowanie pozwala sformułować następujące twierdzenie: Strona 25.

(28) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. TWIERDZENIE 6 Dla każdego układu węzłów siatki można znaleźć ciągłą funkcję f (x) określoną w przedziale [a, b], dla której ciąg wielomianów interpolacyjnych odchyla się znacząco na odcinku [a, b].. Oszacowanie stałej Lebesgue’a przedstawione powyżej (wzory (1.54) oraz (1.55)) odnosi się do przypadku, w którym stosowany jest dowolny zbiór węzłów interpolacyjnych, spełniających warunek (1.14). Dokładność interpolacji (1.55) pogarsza się w przypadku zastosowania konkretnych schematów generowania węzłów w przedziale [a, b]. Najczęściej stosowanym układem punktów interpolacji jest zbiór węzłów równoodległych w przedziale [a, b]. Jest to spowodowane faktem, iż taką siatkę można z łatwością zastosować w numerycznych rozważaniach dowolnych procesów oraz w urządzeniach rejestrujących wyniki eksperymentów. Niestety oszacowanie wartości stałej Lebesgue’a dla tej siatki jest dużo gorsze niż w przypadku dowolnie dobieranych siatek. Turetskii [51] pokazał, że stała Lebesgue’a dla interpolacji w przedziale [-1, 1] jest równa asymptotycznie następującemu wyrażeniu:. SNE (I ) . 2 N 1 ,N , eN ln( N ). (1.56). gdzie: SE oznacza siatkę o węzłach równoodległych w przedziale [a, b]:. i 1   S E   xi : xi  a  (b  a) , i  1, 2, , N  . N 1  . (1.57). Dokładniejsze oszacowania zostały dowiedzione w pracach [50, 52, 53] jednak wyrażenia w nich prezentowane nieznacznie poprawiają (1.56). Szybkość wzrostu wartości stałej Lebesgue’a w przypadku węzłów równoodległych jest znacząco większa niż w przypadku optymalnej siatki. Wybór sitaki o węzłach rozmieszczonych równomiernie w analizowanym przedziale nie gwarantuje wysokiej dokładności interpolacji. Jak to było już jednak wspomniane, siatka ta jest uniwersalna i można ją z powodzeniem stosować podczas analizy rzeczywistych procesów jak i podczas numerycznego rozwiązywania zagadnień inżynierskich. Wartość dużo bliższą optymalnej (1.55) mają stałe Lebesgue'a dla siatek opartych na węzłach pochodzących od wielomianów Czebyszewa [44, 48, 54]. Wielomiany Czebyszewa I rodzaju zdefiniowane są wzorem rekurencyjnym:. Strona 26.

(29) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. T0 ( x)  1,  T1 ( x)  x, T ( x)  2 xT ( x)  T ( x), n  1,2, , N n n1  n1. (1.58). Zachowanie funkcji Lebesgue’a dla węzłów skojarzonych z pierwiastkami lub ekstremami tak zdefiniowanych wielomianów jest dużo bardziej obiecujące. Ograniczenie dolne wartości stałej Lebesgue’a jest na poziomie zbliżonym do wartości dla wielomianu optymalnego. Niech w przedziale [-1, 1] będą dane następujące siatki:    2i  1  ST   xi : xi  cos  , i  1, 2, 3,, N  ,  N   . ST *.    2i  1  cos       N  , i  1, 2, 3,  , N  ,   xi : xi        cos    2 N    .   1  i 1  SCGL   xi : xi  1  cos   , i  1, 2, 3, , N  , 2  N 1   . (1.59a). (1.59b). (1.59c). można pokazać, że dla każdej z wymienionych powyżej siatek stała Lebesgue’a jest ograniczona od dołu przez wartość wzrastającą nieograniczenie wraz z N [44, 55, 56]:. min. S{ST ,ST * ,SCGL }. SN ( I )  2 ln( N )  2 (  ln 8  23 )  O( ln1N ) .. (1.60). Wartość minimalna dla prezentowanych siatek, pochodzących od zer wielomianów Czebyszewa, osiągana jest dla siatki (1.59b). Można zauważyć, że wartość stałej Lebesgue'a jest znacznie mniejsza dla każdej spośród siatek pochodzących od wielomianów Czebyszewa (1.58), niż w przypadku siatki o węzłach rozmieszczonych równomiernie w przedziale (1.57). Węzły tych siatek są znacznie gęściej położone w pobliżu brzegu. Własność ta obniża wartość stałej Lebesgue'a, jednak powoduje równocześnie, że siatki te są dużo mniej uniwersalne w zastosowaniach praktycznych. Dlatego stosowane były tylko do specyficznych zagadnień: przede wszystkim analizy drgań.. 1.3.3 Teoretyczna analiza dokładności przybliżenia pochodnej metodą kwadratur różniczkowych Analiza błędu wynikającego z aproksymacji funkcji i jej pochodnej przy użyciu metody KR została przeprowadzona w pracach doktorskich Shu [28] oraz kontynuowana w rozprawie Chena [57], której promotorem był Shu. Najważniejsze wnioski opracowane przez tych naukowców zostały później zgromadzone przez Shu w monografii [1], na podstawie której w dużej mierze opiera się analiza przeprowadzona w tym paragrafie.. Strona 27.

(30) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. Przeprowadzone przez nich obliczenia wskazują, że dokładność przybliżenia związana jest z wartościami jakie przyjmuje wielomian M (x) dla wybranej siatki. 1.3.3.1 Przybliżenie funkcji wielomianem interpolacyjnym W paragrafie tym zostanie oszacowany błąd aproksymacji funkcji f (x) określonej w przedziale [a, b] wielomianem stopnia (N-1). Zastosowana metoda przybliżenia musi zapewniać funkcji aproksymującej zachowanie warunku interpolacji (1.44). Zadanie takie spełnia wielomian interpolacyjny stopnia (N-1). Stosując interpolację w bazie wielomianów Lagrange'a, na N - węzłowej siatce S, otrzymuje się wielomian interpolacyjny postaci (1.44). Niech E (f ) oznacza błąd interpolacji wielomianem (N-1) stopnia, zdefiniowany następująco: E  f   f  I NS  f  .. (1.61). Rozważany jest przypadek, gdy pochodna N-tego rzędu funkcji f (x) nie jest stała, ale jest ograniczona w przedziale [a, b]. W takiej sytuacji błąd przybliżenia funkcji wielomianem interpolacyjnym, E (f ), może zostać oszacowany przy użyciu wielomianu M (x) w sposób podany poniżej. Niech dana będzie funkcja F (z), zdefiniowana przy pomocy wielomianu interpolacyjnego Lagrange'a (1.44) oraz wielomianu M (x) (1.36): F z   f z   I NS  f z   c M z  .. (1.62). Jasnym jest, że dla tak zdefiniowanej funkcji F (z) w węzłach interpolacji xi (i=1, 2,...,N) występują miejsca zerowe: F xi   0, i  1, 2,, N .. (1.63). Stała c   występująca w funkcji F (z) zdefiniowanej jak (1.62) może być dopasowana dowolnie. Poza punktami xi (i=1, 2, ..., N) różnica f z   I NS  f z  jest różna od zera. Poszukiwana jest taka wartość stałej c  , by: f z   I NS  f z   c M z  .. (1.64). Sytuacja taka ma miejsce dla pewnego x0  (a, b), takiego, że F (x0) = 0. Przypadek taki jest możliwy przy pewnej wartości stałej c. Można stąd wnioskować, że F (z) ma (N+1) pierwiastków w przedziale [a, b]. Stosując N - krotnie twierdzenie Rolle'a do funkcji F (z) można dojść do wniosku, iż F(N) (z) posiada co najmniej jeden pierwiastek znajdujący się w przedziale (a, b). Oznaczając ten pierwiastek przez   (a, b) otrzymuje się:. Strona 28.

(31) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. F  N     0.. (1.65). N S N    i I f z  0 M z   N! , przekształcając (1.65) zapisane jak (1.64) N z N N z otrzymuje się wartość stałej c, przy której F (z) zeruje się w z = x0:. Skoro. f  N    . c N!. (1.66). Zatem błąd aproksymacji wielomianem interpolacyjnym określa zależność: E f  . f  N    M x  . N!. (1.67). W powyższym równaniu zmienna  jest funkcją zależną od zmiennej x. 1.3.3.2 Związek KR z pochodną wielomianu interpolacyjnego Zgodnie z definicją kwadratury różniczkowej (1.15) przybliżenie m-tej pochodnej funkcji f (x) w punkcie xi  S (1.14) odbywa się według wzoru: N. *. N. f m  xi    i,mj  f j   f j lkm  xi , m  1, 2, , N  1, i  1, 2, , N , j 1. (1.68). j 1. gdzie: równość (*) jest spełniona na mocy (1.34). Pochodna wielomianu interpolacyjnego Lagrange'a (1.44) jest równa:. m m p N x   m x m x. N. N.  f l x   f l   x, m. k 1. k k. k 1. k k. m  1, 2, , N  1.. (1.69). Zestawiając ze sobą wzory (1.68) oraz (1.69) można sformułować wniosek, iż: Przybliżenie m-tej pochodnej funkcji f w punkcie xi  S (1.14) przy pomocy KR jest równe wartości odpowiedniej pochodnej wielomianu interpolacyjnego Lagrange'a w punkcie xi zbudowanego na siatce S. 1.3.3.3 Przybliżenie dowolnej pochodnej funkcji kwadraturą różniczkową Podstawy do analiz przeprowadzonych w tym paragrafie opierają się, podobnie jak w przypadku paragrafu §1.3.3.1, na pracach [1, 28, 57]. W związku z wynikami uzyskanymi w poprzednim paragrafie, analizy metody KR przeprowadzone tylko w węzłach interpolacji prowadzić mogą do zaniedbań istotnych wartości. Skoro współczynniki wagowe metody KR związane są z pochodną wielomianu interpolacyjnego, istotne znaczenie ma jak szybko wartości tego wielomianu oddalają się od wartości aproksymowanej funkcji. Niech ED( m)  f  oznacza błąd przybliżenia m-tej pochodnej funkcji f pochodną wielomianu interpolacyjnego: Strona 29.

(32) Zastosowanie metody kwadratur różniczkowych w komputerowej symulacji przewodzenia ciepła i krzepnięcia odlewów. EDm   f  . . .  m f  m I NS  f  m   f  I NS  f  . add x m x m x m. . . (1.70). W powyższym wzorze m = 1, 2, ..., N-1, add - oznacza, iż w miejscu tym korzysta się z własności addytywności operatora różniczkowania. Zgodnie z wnioskiem płynącym z poprzedniego paragrafu (§1.3.3.2) błąd ten w węzłach siatki jest równy błędowi przybliżenia pochodnej metodą KR. Na podstawie (1.67): 1 m f  N    M x  . m N! x. . EDm   f  . . (1.71). Wyrażenie to pozwala wyznaczyć błędy przybliżenia kolejnych pochodnych funkcji f pochodnymi wielomianu interpolacyjnego: ED  f   1. ED2   f  . f  N 1  M x  1  f  N   M 1 x  , N!. .  . 1  N 2   M x   1 2  2 f  N 1  M 1 x  1  f N!  f  N 1  M x  2   f  N   M 2  x  ,. . E D3  f  . .  .  M x   f  M x    f  N 2   M 2  x  1  f  N 2   M 2  x  2   f  N   M 3 x . 2f. 2 . 1.  N 1. (1.73).  . 1  N 3   M x   1 3  3 f  N 2   M 1 x   1  f N!  3 f  N 1  M x  1 2   2 f  N 1  M 1 x  2    N 1. (1.72). 3 . (1.74). W analogiczny sposób można wyznaczyć wartość błędu dla pochodnych kolejnych rzędów. Powyższa analiza pozwala stwierdzić, że błąd przybliżenia m-tej pochodnej funkcji f m-tą pochodną wielomianu interpolacyjnego zależy od przybliżanej funkcji f, wartości x oraz wartości wielomianu M (x) oraz jego m kolejnych pochodnych. Wartość wielomianu M (x) jest nieograniczona, gdy wzrasta liczba węzłów siatki (§1.3.2.1). Od pewnej liczby punktów siatki wzrost jakości przybliżenia funkcji o skomplikowanych kształtach spowodowany wzrostem dokładności odwzorowania jej kształtu parami punktów (xi, fi), będzie zniwelowany przez wzrost wartości stałej Lebesgue'a dla interpolacji. Ponadto można spodziewać się, iż wartość pochodnych wielomianu M (x) będzie jeszcze większa. Nasuwa się zatem wniosek, że począwszy od siatek zawierających pewną, znaczną liczbę węzłów interpolacji (dyskretyzacji) dalsze zagęszczanie siatki spowoduje wzrost błędu przybliżenia.. Strona 30.