1. Rozwiązać równanie:
z warunkiem początkowym y(1) = y'(1) = 0 w przedziale x ∈ (1, 5) najpierw metodą Eulera z krokiem h = 0,1, potem analitycznie. Rozwiązania przedstawić na jednym wykresie.

Napisać równanie wielomianu Lagrange'a i Newtona, interpolującego funkcję y(x) w 5 węzłach rozłożonych równomiernie.

Zgodnie z teorią podaną na wykładzie, dla równań różniczkowych typu

(1)

stosujemy następujące przejścia: y → y⁽⁰⁾,
,
, …,
,
. Wtedy równanie różniczkowe przyjmuje postać:

(1a)

a poszczególnych zmiennych szukamy z układu równań:

(2)

W programie zamiast y⁽⁰⁾, y⁽¹⁾, y⁽²⁾, … piszemy y0, y1, y2, …

Jeśli równanie (1) rozwiązujemy metodą Eulera, stosujemy algorytm

(3)

a więc układ (2) będzie miał następującą postać numeryczną:

(4)

Równanie w zadaniu (gdzie n = 2) przekształci się następująco:

⇒

⇒
(5)

Program zaczynamy od wypisania wektora argumentów (przedział (1, 5), skok h = 0,1):

In[1]:= x0 = Table[i,{i,1,5,.1}]

Out[1]:= {1.,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.,2.1,2.2,2.3,

2.4,2.5,2.6,2.7,2.8,2.9,3.,3.1,3.2,3.3,3.4,3.5,3.6,3.7,

3.8,3.9,4.,4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5.}

Następnie - szablony dla wektorów y0 i y1 (takiego samego rozmiaru, co wektor x0):

In[2]:= y0 = Table[0,{i,1,Length[x0]}]

Out[2]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0}

In[3]:= y1 = Table[0,{i,1,Length[x0]}]

Out[3]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0}

Warunki początkowe: y(1) = y'(1) = 0 ⇒ y⁽⁰⁾|_x=1 = 0, y⁽¹⁾|_x=1 = 0. Oznacza to, że trzeba znaleźć element wektora x0 równy 1 i elementowi wektora y0 i y1 o takim samym numerze przyporządkować wartość 0. W wektorze x0 jest to element 1-szy (tak jest najczęściej, ale nie zawsze), a więc ma być
,
:

In[4]:= y0[[1]] = 0

Out[4]:= 0

In[5]:= y1[[1]] = 0

Out[5]:= 0

Uwaga: Zapis y0[[1]] = 0 nie oznacza, że y dla x = 1 wynosi 0 (bo tak można by się zasugerować po przeczytaniu treści zadania), ale to, że pierwsza wartość y wynosi 0.

Teraz można pisać pętlę (taką, jak na zajęciach). Zawieramy w niej równania z układu (5).

In[6]:= Do[y0[[i+1]]=y0[[i]]+.1 y1[[i]]; y1[[i+1]]=y1[[i]]+

.1
,{i,1,Length[x0]-1}]

(gdybyśmy napisali Length[x0], indeksowanie y0[[i+1]] - zgodne z układem (5) − wymusiłoby na pętli policzenie jednego wyrazu więcej, niż wynosi długość wektora x0)

Rezultat otrzymamy po napisaniu nazwy wektora, którego szukamy - y0:

In[7]:= y0

Out[7]:= {0,0,0.005,0.015,0.0300104,0.0500497,0.0751425,

0.105319,0.140612,0.181061,0.226705,0.277587,0.333752,

0.395248,0.462124,0.534428,0.612215,0.695535,0.784445,

0.878998,0.979252,1.08526,1.19709,1.31479,1.43843,

1.56807,1.70376,1.84557,1.99356,2.14781,2.30836,2.47529,

2.64866,2.82854,3.015,3.20811,3.40792,3.61452,3.82797,

4.04835,4.27571}

Jest to rozwiązanie numeryczne. Rozwiązanie analityczne dostajemy, używając funkcji DSolve, pamiętając o tym, aby:

− pochodną zapisywać jako „primy” (np.
→ y''),

− pisać podwójny znak równości,

− szukaną funkcję pisać razem z argumentem (nie samo y, ale y[x]).

Warunki początkowe zapisujemy razem z równaniem w nawiasie klamrowym. Potem piszemy, jakiej funkcji szukamy (y[x]) i jaki jest jej argument (x).

In[8]:= DSolve[{y''[x]-.25
.5,y[1]0,y'[1]0},y[x],x]

Out[8]:= {{y[x]→ -0.750322
BesselI[1.,
]+0.271495
BesselK

[1.,
]-0.5x^5/2BesselK[1.,
]Hypergeometric0F1Regulari

zed[3.,0.25 x]+8.BesselI[1.,
]MeijerG[{{1.5},{}},

{{1.5,2.5},{0.5}},0.5
,0.5]}}

Rozwiązanie jest nie do rozszyfrowania, ale nie należy się tym przejmować. Po narysowaniu wykresu otrzymanej funkcji wraz z wykresem rozwiązania numerycznego będziemy mieli odpowiedź, w jakim stopniu one sobie odpowiadają (oczywiście, powinny się jak najlepiej pokrywać). Aby sprawniej operować rozwiązaniem analitycznym, wczytamy je w pamięć Mathematiki jako funkcję:

In[9]:= f[x_] -0.750322
BesselI[1.,
]+0.271495
BesselK

[1.,
]-0.5x^5/2BesselK[1.,
]Hypergeometric0F1Regulari

zed[3.,0.25 x]+8.BesselI[1.,
]MeijerG[{{1.5},{}},

{{1.5,2.5},{0.5}},0.5
,0.5]

(nie ma odpowiedzi, bo przed „=” stoi dwukropek). Wykresy:

In[10]:= Plot[f[x],{x,1,5},Epilog→{Table[Point[{x0[[i]],y0[[i]]}],

{i,1,Length[x0]}]}]

Out[10]:=

Teraz przybliżymy funkcję f(x) (wczytaną w instrukcji 9) wielomianem Lagrange'a i Newtona w 5 węzłach - a więc wielomian przybliżający będzie 4. stopnia. Węzły:

In[11]:= xr = {1,2,3,4,5}

Out[11]:= {1,2,3,4,5}

Wielomian Lagrange'a:
(u nas N = 4). Warunek m ≠ k oznacza, że w iloczynie nie będzie ułamka
, co można zapisać funkcją warunkową: jeśli m ≠ k, to ułamek wystąpi w postaci podanej we wzorze, jeśli m = k - będzie równy 1 (bo 1 nie zmienia wartości iloczynu). Czyli:

(6)

Trzeba też jeszcze pamiętać, że nie możemy indeksować od 0, a więc numerowanie k i m zaczynamy od 1, dlatego kończymy nie na 4, tylko na 5. Polecenie napisania wielomianu oraz jego rezultat są następujące:

In[12]:=

Out[12]:= 6.39718×10^-7(2-x)(3-x)(4-x)(5-x)+0.0422265(3-x)(4-x)

(5-x)(-1+x)+0.26001(4-x)(5-x)(-2+x)(-1+x)+0.402355

(5-x)(-3+x)(-2+x)(-1+x)+0.184986(-4+x)(-3+x)(-2+x)(-1+x)

W postaci sumy:

In[13]:= Expand[WL]

Out[13]:= 0.235907-0.468301x+0.227135x²+0.00485823x³+0.000415378x⁴

Wielomian Newtona:

(7a), gdzie
(7b)

(u nas N = 4). Indeksy przy b są różne: indeks n we wzorze (7a) oznacza, że wyrazów b jest tyle, ile składników w sumie w (7a) i każdemu z nich przyporządkowany jest inny wyraz b, z kolei indeks m we wzorze (7b) oznacza, że składników w sumie w (7b) jest tyle, ile wynosi numer wyrazu b. Nie można tych oznaczeń pomylić ani przyjąć jako jedno - jak widać, w (7b) jest również indeks n i oznacza on dokładnie to, co w (7a), czyli numer kolejnego składnika w sumie w (7a).

Warunek v ≠ n oznacza - jak poprzednio - że w iloczynie nie będzie czynnika x_n - x_v, więc

(7c)

Podobnie, jak wcześniej, indeksowanie zaczynamy nie od 0, ale od 1 - dlatego we wz. (7a) n zmienia się od 1 do 5 (nie od 0 do 4), a v - od 2 (nie od 1). Z kolei we wz. (7c) n i v zmienia się od 1, nie od 0.

Najpierw trzeba obliczyć wyrazy b. Zależnie od wartości indeksu, mamy tu różną liczbę składników (gdy rośnie indeks, dochodzi jeden składnik), dlatego do obliczenia stosujemy pętlę (Do), w której numer kroku wyznacza liczbę składników w sumie. Do tego będzie potrzebny szablon na b - będzie on zawierał 5 wyrazów, bo we wz. (7a) n zmienia się od 1 do 5:

In[14]:= b = Table[0,{i,1,5}]

Out[14]:= {0,0,0,0,0}

In[15]:= Do[
,{m,1,5}]

In[16]:= b

Out[16]:= {0.0000153532,0.253343,0.266669,0.00901201,0.000415378}

Teraz liczymy wielomian:

In[17]:=

Out[17]:= 0.0000153532+0.253343(-1+x)+0.266669(-2+x)(-1+x)+

0.00901201(-3+x)(-2+x)(-1+x)+0.000415378(-4+x)(-3+x)

(-2+x)(-1+x)

W postaci sumy - wzór jest ten sam, co dla WL:

In[18]:= Expand[WN]

Out[18]:= 0.235907-0.468301x+0.227135x²+0.00485823x³+0.000415378x⁴

Wszystkie sumy można zapisać również jako Sum, a iloczyny jako Product. Czyli np. wiersz 17 może mieć też postać:

In[17]:= WN = Sum[b[[n]] Product[x-xr[[v-1]],{v,2,n}],{n,1,5}]

Instrukcja Plot pokaże, że wielomiany bardzo dobrze przybliżają funkcję f(x):

In[19]:= Plot[{f[x],WL},{x,1,5},PlotStyle→{Thickness[.004],

Thickness[.008]}]

Out[19]:=

Wykresy się pokrywają (tu - f(x) i WL).

2. Rozwiązać równanie:
z warunkiem początkowym y(−1) = 1 w przedziale x ∈ (−1, 3) najpierw metodą Heuna z krokiem h = 0,2, potem analitycznie. Rozwiązania przedstawić na jednym wykresie.

Napisać równania splajnu 3. stopnia interpolującego funkcję y(x) w 4 węzłach: −1, 0, 1, 2.4, oraz splajnu 2. stopnia interpolującego funkcję y(x) w 5 węzłach: −0.5, 0.6, 1.5, 2.3, 3. Narysować każdy ze splajnów wraz z funkcją y(x).

Wg teorii podanej na wykładzie, algorytm metody Heuna dla równania różniczkowego
jest następujący:

− pierwszy predyktor:

− pierwsze przybliżenie:

− korektor:

− drugi predyktor:

− właściwe przybliżenie - szukane rozwiązanie:

Po sprowadzeniu równania do postaci odpowiedniej dla obliczeń numerycznych (pochodna po jednej stronie, reszta - po drugiej stronie znaku równości) mamy:
, a więc
. Zatem:

− pierwszy predyktor:

− pierwsze przybliżenie:

− korektor:
(8)

− drugi predyktor:

− właściwe przybliżenie - szukane rozwiązanie:

Program zaczynamy od wypisania wektora argumentów (przedział (−1, 3), skok h = 0,2):

In[1]:= x0 = Table[i,{i,−1,3,.2}]

Out[1]:= {-1.,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.,1.2,1.4,

1.6,1.8,2.,2.2,2.4,2.6,2.8,3.}

Następnie - szablony dla poszczególnych wielkości z algorytmu, czyli p1, y0, k, p2, y1 (tę ostatnią oznaczamy y1, a nie y, dlatego, że przez y będzie oznaczona funkcja w poleceniu DSolve):

In[2]:= y1=Table[0,{i,1,Length[x0]}]

Out[2]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

In[3]:= p1=Table[0,{i,1,Length[x0]}]

Out[3]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

In[4]:= y10=Table[0,{i,1,Length[x0]}]

Out[4]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

In[5]:= k=Table[0,{i,1,Length[x0]}]

Out[5]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

In[6]:= p2=Table[0,{i,1,Length[x0]}]

Out[6]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

Warunek początkowy: szukamy elementu wektora x0 równego −1 i elementowi wektora y1 o takim samym numerze przyporządkowujemy wartość 1. W wektorze x0 jest to element 1-szy (tak jest najczęściej, ale nie zawsze), a więc ma być y₁ = 1:

In[7]:= y1[[1]] = 1

Out[7]:= 1

Teraz piszemy pętlę (taką, jak na zajęciach). Zawieramy w niej równania z algorytmu (8):

In[8]:= Do[p1[[i]]=Sin[3x0[[i]]]^.5x0[[i]]-y1[[i]]; y10[[i]]=

y1[[i]]+.2 p1[[i]]; k[[i]]=Sin[3x0[[i+1]]]^{.5 x0[[i+1]]}-

y10[[i]]; p2[[i]]=
; y1[[i+1]]=y1[[i]]+

.2 p2[[i]],{i,1,Length[x0]-1}]

Rezultat otrzymamy po napisaniu nazwy wektora, którego szukamy - y1:

In[9]:= y1

Out[9]:= {1,0.767875,0.521291,0.293434,0.128478,0.0644791,

0.115275,0.258288,0.434323,0.562077,0.564783,0.401103,

0.088885,-0.289226,-0.604595,-0.723776,-0.560667,

-0.121366,0.477681,1.03137,1.30767}

Jest to rozwiązanie numeryczne. Rozwiązanie analityczne dostajemy, używając funkcji DSolve:

In[10]:= DSolve[{y'[x]+y[x] ^.5xSin[3x],y[-1]1},y[x],x]

Out[10]:= {{y[x]→

}}

Aby sprawniej posługiwać się otrzymanym rozwiązaniem, wczytamy je w pamięć jako funkcję f(x):

In[11]:= f[x_]=Simplify[

]

Out[11]:= e^-1-x +

(nie ma dwukropka po znaku równości - jest odpowiedź).

Wykresy:

In[12]:= Plot[f[x],{x,-1,3},Epilog→{Table[Point[{x0[[i]],

y1[[i]]}],{i,1,Length[x0]}]}]

Out[12]:=

Teraz przybliżymy funkcję f(x) (wczytaną w instrukcji 11) splajnami. Pierwszy to splajn 3. stopnia, poprowadzony w 4 węzłach: x₀ = −1, x₁ = 0, x₂ = 1, x₃ = 2,4. Wszystkie węzły dobieramy tak, żeby przybliżenie było jak najlepsze (tu również węzły nie zostały dobrane przypadkowo, ale w wyniku kilku prób). Dlatego też węzeł pierwszy i ostatni nie muszą pokrywać się z granicą przedziału.

In[13]:= x1 = {-1,0,1,2.4}

Out[13]:= {-1,0,1,2.4}

Splajn 3. stopnia, poprowadzony w 4 węzłach jest opisany następującymi równaniami:

f₀ = a₀(x − x₀)³ + b₀(x − x₀)² + c₀(x − x₀) + d₀

f₁ = a₁(x − x₁)³ + b₁(x − x₁)² + c₁(x − x₁) + d₁ (9)

f₂ = a₂(x − x₂)³ + b₂(x − x₂)² + c₂(x − x₂) + d₂

Na wykładzie były podane macierze układu równań, służącego do wyznaczenia współczynni-ków splajnu 3. stopnia poprowadzonego w 5 węzłach. Macierz A pokazano na następnej stronie (trzeba pamiętać, że ponieważ nie ma indeksu zerowego, numerujemy węzły od 1 - dlatego też numeracja kończy się na 5). Z tej macierzy wykreślamy odpowiednie kolumny i wiersze - zależnie od tego, jaki mamy stopień wielomianu i przez ile węzłów jest on poprowadzony (można by też dodawać kolumny i wiersze - gdyby było więcej węzłów niż 5 albo wyższy stopień niż 3 - ale tym się nie zajmujemy). Ponieważ liczymy splajn 3. stopnia, poprowadzony w 4 węzłach, wykreślamy te kolumny, które zawierają węzeł 5-ty, czyli 4., 8. i 11. kolumnę oraz odpowiadające im wiersze (4., 8., 11.). Ze względu na brak miejsca, poniżej znajdują się tylko wyniki obliczeń.

In[14]:= A = ...

Out[14]:= {{1,0,0,1,0,0,0,0},{0,1,0,0,1,0,1,0},{0,0,2.744,0,0,

1.96,0,1.4},{3,0,0,2,0,0,-1,0},{0,3,0,0,2,0,1,-1},

{0,0,5.88,0,0,2.8,0,1},{6,0,0,2,-2,0,0,0},{0,6,0,0,2,

-2,0,0}}

Postać macierzowa (nie musi być, ale pozwala na sprawdzenie, czy macierz jest kwadratowa):

In[15]:= MatrixForm[A]

Out[15]//MatrixForm=

Macierz w przekształcamy następująco:

− wykreślamy wiersz 4., 7. i 11.,

− w wierszu 8. zamiast węzła 5 wpisujemy węzeł 4.

Wynik:

In[16]:= w = ...

Out[16]:= {{0.132099},{0.534986},{-0.740351},{1.08559},{0.},

{2.79393},{0.},{0.}}

Teraz obliczamy współczynniki równań (9) - z równania macierzowego A wsp = w (gdzie wsp - szukany wektor współczynników). Wektor wsp = A⁻¹w zawiera tylko współczynniki a₀, a₁, a₂, b₀, b₁, b₂, c₁, c₂. Pozostałe są równe:

c₀ = f `(x₀) = D[f[x],x]/.x→x1[[1]] d₀ = f(x₀) = f[x1[[1]]]

d₁ = f(x₁) = f[x1[[2]]] d₂ = f(x₂) = f[x1[[3]]]

In[17]:= wsp=Inverse[A].w//N

Out[17]:= {{0.895428},{-1.462},{1.7273},{-0.763329},{1.92296},

{-2.46305},{0.0740323},{-0.466063}}

Równania (9) przyjmą postać:

f₀ = f(x₀) + f `(x₀) (x − x₀) + b₀(x − x₀)² + a₀(x − x₀)³

f₁ = f(x₁) + c₁(x − x₁) + b₁(x − x₁)² + a₁(x − x₁)³ (9a)

f₂ = f(x₂) + c₂(x − x₂) + b₂(x − x₂)² + a₂(x − x₂)³

Kolejność jest zmieniona po to, aby wygodniej wyrazić zmienność potęgi różnicy x − x_k. Równania (9a) można przedstawić w postaci iloczynu skalarnego:

f₀ = [(x − x₀)⁰, (x − x₀)¹, (x − x₀)², (x − x₀)³]*[f(x₀), f `(x₀), b₀, a₀]

f₁ = [(x − x₁)⁰, (x − x₁)¹, (x − x₁)², (x − x₁)³]*[f(x₁), c₁, b₁, a₁] (9b)

f₂ = [(x − x₂)⁰, (x − x₂)¹, (x − x₂)², (x − x₂)³]*[f(x₂), c₂, b₂, a₂]

Ponieważ teraz potęga rośnie, pierwsze wektory można zapisać w postaci tablicy:

Table[(x-x1[[1]])^k,{k,0,3}], Table[(x-x1[[2]])^k,{k,0,3}],...

Aby obliczyć funkcje f₀ ÷ f₂, do równań (9b) wstawiamy współczynniki otrzymane w wektorze wsp - w kolejności pokazanej wcześniej (trzeba pamiętać, że jeśli wektor jest kolumnowy, należy go potraktować jako macierz, a więc wyciągać element, podając numer zarówno wiersza, jak i kolumny - w tym przypadku 1): a₀ - wsp[[1,1]],

a₁ − wsp[[2,1]], …, c₁ − wsp[[7,1]], c₂ − wsp[[8,1]]:

In[18]:= f0=Table[(x-x1[[1]])^k,{k,0,3}].{f[x1[[1]]],

D[f[x],x]/.x→x1[[1]],wsp[[4,1]],wsp[[1,1]]}//N

Out[18]:= 1.-1.08559(1.+x)-0.763329(1.+x)²+0.895428(1.+x)³

In[19]:= f1=Table[(x-x1[[2]])^k,{k,0,3}].{f[x1[[2]]],wsp[[7,1]],

wsp[[5,1]],wsp[[2,1]]}//N

Out[19]:= 0.0465053+0.0740323x+1.92296x²-1.462x³

In[20]:= f2=Table[(x-x1[[3]])^k,{k,0,3}].{f[x1[[3]]],wsp[[8,1]],

wsp[[6,1]],wsp[[3,1]]}//N

Out[20]:= 0.581491-0.466063(-1.+x)-2.46305(-1.+x)²+1.7273(-1.+x)³

Łączymy otrzymane funkcje w całość za pomocą polecenia If:

In[21]:= int[x_]=If[x<xa[[2]],f0,If[x<xa[[3]],f1,f2]]

Out[21]:= If[x<0,f0,If[x<xa[[3]],f1,f2]]

Wykres (linia gruba - interpolant, cienka - funkcja interpolowana):

In[22]:= Plot[{f[x],int[x]},{x,-1,3},GridLines→{x1,{0}},

PlotStyle→{Thickness[.003],Thickness[.006]}]

Out[22]:=

Teraz obliczymy splajn 2. stopnia, poprowadzony w 5 węzłach: x₀ = −0,5, x₁ = 0,6, x₂ = 1,5, x₃ = 2,3, x₄ = 3:

In[23]:= x2 = {-0.5,0.6,1.5,2.3,3}

Out[23]:= {-0.5,0.6,1.5,2.3,3}

Splajn 2. stopnia, poprowadzony w 5 węzłach jest opisany następującymi równaniami:

f₀ = a₀(x − x₀)² + b₀(x − x₀) + c₀

f₁ = a₁(x − x₁)² + b₁(x − x₁) + c₁ (10)

f₂ = a₂(x − x₂)² + b₂(x − x₂) + c₂

f₃ = a₃(x − x₃)² + b₃(x − x₃) + c₃

Z macierzy A wykreślamy cztery pierwsze kolumny i 4 ostatnie wiersze.

In[24]:= A = ...

Out[24]:= {{1.21,0,0,0,0,0,0},{0,0.81,0,0,0.9,0,0},{0,0,0.64,0,

0,0.8,0},{0,0,0,0.49,0,0,0.7},{2.2,0,0,0,-1,0,0},

{0,1.8,0,0,1,-1,0},{0,0,1.6,0,0,1,-1}}

Postać macierzowa:

In[25]:= MatrixForm[A]

Out[25]//MatrixForm=

Macierz w przekształcamy, wykreślając 4 ostatnie wyrazy:

Wynik:

In[26]:= w = ...

Out[26]:= {{1.32327},{-0.515976},{-0.324972},{1.76278},{1.17491},

{0.},{0.}}

Obliczamy współczynniki równań (10) z równania macierzowego A wsp = w ⇒ wsp = A⁻¹w. Wektor wsp zawiera tylko współczynniki a₀, a₁, a₂, a₃, b₁, b₂, b₃. Pozostałe są równe:

b₀ = f `(x₀) = D[f[x],x]/.x→x2[[1]] c₀ = f(x₀) = f[x2[[1]]]

c₁ = f(x₁) = f[x2[[2]]] c₂ = f(x₂) = f[x2[[3]]]

c₃ = f(x₃) = f[x2[[4]]]

In[27]:= wsp=Inverse[A].w//N

Out[27]:= {{1.09361},{-2.00483},{2.46429},{1.36148},{1.23104},

{-2.37765},{1.56522}}

Równania (10) przyjmą postać:

f₀ = f(x₀) + f `(x₀) (x − x₀) + a₀(x − x₀)²

f₁ = f(x₁) + b₁(x − x₁) + a₁(x − x₁)² (10a)

f₂ = f(x₂) + b₂(x − x₂) + a₂(x − x₂)²

f₃ = f(x₃) + b₃(x − x₃) + a₃(x − x₃)²

Równania (10a) można przedstawić w postaci iloczynu skalarnego:

f₀ = [(x − x₀)⁰, (x − x₀)¹, (x − x₀)²]*[f(x₀), f `(x₀), a₀]

f₁ = [(x − x₁)⁰, (x − x₁)¹, (x − x₁)²]*[f(x₁), b₁, a₁] (10b)

f₂ = [(x − x₂)⁰, (x − x₂)¹, (x − x₂)²]*[f(x₂), b₂, a₂]

f₃ = [(x − x₃)⁰, (x − x₃)¹, (x − x₃)²]*[f(x₃), b₃, a₃]

Pierwsze wektory można zapisać w postaci tablicy:

Table[(x-x2[[1]])^k,{k,0,2}], Table[(x-x2[[2]])^k,{k,0,2}],...

Aby obliczyć funkcje f₀ ÷ f₃, do równań (10b) wstawiamy współczynniki otrzymane w wektorze wsp: a₀ - wsp[[1,1]], a₁ − wsp[[2,1]], …, b₂ − wsp[[6,1]],

b₃ − wsp[[7,1]]:

In[28]:= f0=Table[(x-x2[[1]])^k,{k,0,2}].{f[x2[[1]]],

D[f[x],x]/.x→x2[[1]], wsp[[1,1]]}//N

Out[28]:= 0.398063-1.17491(0.5+x)+1.09361(0.5+x)²

In[29]:= f1=Table[(x-x2[[2]])^k,{k,0,2}].{f[x2[[2]]],

wsp[[5,1]],wsp[[2,1]]}//N

Out[29]:= 0.428931+1.23104(-0.6+x)-2.00483(-0.6+x)²

In[30]:= f2=Table[(x-x2[[3]])^k,{k,0,2}].{f[x2[[3]]],

wsp[[6,1]],wsp[[3,1]]}//N

Out[30]:= -0.0870448-2.37765(-1.5+x)+2.46429(-1.5+x)²

In[31]:= f3=Table[(x-x2[[4]])^k,{k,0,2}].{f[x2[[4]]],

wsp[[7,1]],wsp[[4,1]]}//N

Out[31]:= -0.412016+1.56522(-2.3+x)+1.36148(-2.3+x)²

Łączymy otrzymane funkcje w całość za pomocą polecenia If:

In[32]:= int[x_]=If[x<xa[[2]],f0,If[x<xa[[3]],f1,If[x<xa[[4]],

f2,f3]]]

Out[32]:= If[x<0,f0,If[x<xa[[3]],f1,If[x<xa[[4]],f2,f3]]]

Wykres (linia gruba - interpolant, cienka - funkcja interpolowana):

In[33]:= Plot[{f[x],int[x]},{x,-1,3},GridLines→{x2,{0}},

PlotStyle→{Thickness[.003],Thickness[.006]}]

Out[33]:=

3. Rozwiązać równanie:
z warunkiem początkowym y(0) = 2, y'(0) = 0, y''(0) = 1, y'''(0) = 0, w przedziale x ∈ (0,3) najpierw metodą Eulera z krokiem h = 0,01, potem analitycznie. Rozwiązania przedstawić na jednym wykresie.

Napisać równania splajnu 3. stopnia interpolującego funkcję y(x) w 4 węzłach oraz splajnu 2. stopnia interpolującego funkcję y(x) w 5 węzłach - dla obu splajnów węzły rozłożone równomiernie. Narysować każdy ze splajnów wraz z funkcją y(x).

Równanie w zadaniu (gdzie n = 4) przekształci się następująco:

⇒

(poszczególne zmienne y⁽ⁱ⁾ otrzymują numery takie, jak rząd pochodnej funkcji szukanej, najwyższa pochodna funkcji szukanej zmienia się w pierwszą pochodną zmiennej y⁽ⁱ⁾ o numerze i = n − 1)

⇒
(11)

Dane wejściowe (x0, y0, y1, y2, y3 − wektory te mają 301 elementów, dlatego większość została pominięta):

In[1]:= x0 = Table[i,{i,0,3,.01}]

Out[1]:= {0.,0.01,0.02,0.03,...,2.98,2.99,3.}

In[2]:= y0 = Table[0,{i,1,Length[x0]}]

Out[2]:= {0,0,0,...,0,0,0}

In[3]:= y1 = Table[0,{i,1,Length[x0]}]

Out[3]:= {0,0,0,...,0,0,0}

In[4]:= y2 = Table[0,{i,1,Length[x0]}]

Out[4]:= {0,0,0,...,0,0,0}

In[5]:= y3 = Table[0,{i,1,Length[x0]}]

Out[5]:= {0,0,0,...,0,0,0}

Warunki początkowe:
,
,
,
:

In[6]:= y0[[1]] = 2

Out[6]:= 2

In[7]:= y1[[1]] = 0

Out[7]:= 0

In[8]:= y2[[1]] = 1

Out[8]:= 1

In[9]:= y3[[1]] = 0

Out[9]:= 0

Pętla:

In[10]:= Do[y0[[i+1]]=y0[[i]]+.01 y1[[i]]; y1[[i+1]]=y1[[i]]+

.01 y1[[i]]; y2[[i+1]]=y2[[i]]+.01 y1[[i]]; y3[[i+1]]=

y3[[i]]+.01(2 x0[[i]] y1[[i]]-y0[[i]]),

{i,1,Length[x0]-1}]

Wynik:

In[11]:= y0

Out[11]:= {2,2,2.0001,2.0003,2.0006,2.001,2.0015,2.0021,2.0028,

...,1.20844,1.16471,1.12012,1.07467,1.02833,0.981093}

Rozwiązanie analityczne:

In[12]:= DSolve[y''''[x]-2x y'[x]+y[x]0,y[0]2,y'[0]0,

y''[0]1,y'''[0]0},y[x],x]

Out[12]:= {{y[x]→
(4 HypergeometricPFQ
,
+

x² HypergeometricPFQ
,
)}}

W celu sprawniejszego operowania rozwiązaniem, wprowadzamy funkcję f[x_]:

In[13]:= f[x]
(4 HypergeometricPFQ
,
+

x² HypergeometricPFQ
,
)

Wykres:

In[13]:= Plot[f[x],{x,0,3},PlotRange→{f[3]//N,Max[y0]},Epilog→

{Table[Point[{x0[[i]],y0[[i]]}],{i,1,Length[x0]}]}]

Out[13]:=

Splajn 3. stopnia w 4 węzłach rozłożonych równomiernie (funkcja f(x) wczytana w instrukcji 13):

− węzły:

In[14]:= x1 = Table[i,{i,0,3,1}

Out[14]:= {0,1,2,3}

− macierz A:

In[15]:= A = ...

Out[15]:= {{1,0,0,1,0,0,0,0},{0,1,0,0,1,0,1,0},{0,0,1,0,0,

1,0,1},{3,0,0,2,0,0,-1,0},{0,3,0,0,2,0,1,-1},

{0,0,3,0,0,2,0,1},{6,0,0,2,-2,0,0,0},{0,6,0,0,2,

-2,0,0}}

In[16]:= MatrixForm[A]

Out[16]//MatrixForm=

− macierz w:

In[17]:= w = ...

Out[17]:= {0.420495},{0.431569},{-2.01341},{0.},{0.},

{-5.02965},{0.}}

− współczynniki równań (9) - wektor wsp:

In[18]:= wsp=Inverse[A].w//N

Out[18]:= {{-0.178281},{-0.295072},{-1.09748},{0.598776},

{0.0639326}, {-0.821282},{0.662708},{-0.0946412}}

− równania (9):

In[19]:= f0=Table[(x-x1[[1]])^k,{k,0,3}].{f[x1[[1]]],

D[f[x],x]/.x→x1[[1]],wsp[[4,1]],wsp[[1,1]]}//N

Out[19]:= 2.+0.598776x²-0.178281x³

In[20]:= f1=Table[(x-x1[[2]])^k,{k,0,3}].{f[x1[[2]]],wsp[[7,1]],

wsp[[5,1]],wsp[[2,1]]}//N

Out[20]:= 2.42049+0.662708(-1.+x)+0.0639326(-1.+x)²-0.295072

(-1.+x)³

In[21]:= f2=Table[(x-x1[[3]])^k,{k,0,3}].{f[x1[[3]]],wsp[[8,1]],

wsp[[6,1]],wsp[[3,1]]}//N

Out[21]:= 2.85206-0.0946412(-2.+x)-0.821282(-2.+x)²-1.09748

(-2.+x)³

− połączenie otrzymanych funkcji w całość:

In[22]:= int[x_]=If[x<xa[[2]],f0,If[x<xa[[3]],f1,f2]]

Out[22]:= If[x<0,f0,If[x<xa[[3]],f1,f2]]

− wykres (linia gruba - interpolant, cienka - funkcja interpolowana):

In[23]:= Plot[{f[x],int[x]},{x,0,3},GridLines→{x1,{0}},

PlotStyle→{Thickness[.003],Thickness[.006]}]

Out[23]:=

Splajn 2. stopnia w 5 węzłach rozłożonych równomiernie:

− węzły:

In[24]:= x1 = Table[i,{i,0,3,0.75}

Out[24]:= {0,0.75,1.5,2.25,3}

− macierz A:

In[25]:= A = ...

Out[25]:= {{0.5625,0,0,0,0,0,0},{0,0.5625,0,0,0.75,0,0},{0,0,

0.5625,0,0,0.75,0},{0,0,0,0.5625,0,0,0.75},{1.5,0,0,

0,-1,0,0},{0,1.5,0,0,1,-1,0},{0,0,1.5,0,0,1,-1}}

In[26]:= MatrixForm[A]

Out[26]//MatrixForm=

− macierz w:

In[27]:= w = ...

Out[27]:= {{0.25559},{0.486566},{-0.0107309},{-1.89277},{0.},

{0.},{0.}}

− współczynniki równań (10) − wektor wsp:

In[28]:= wsp=Inverse[A].w//N

Out[28]:= {{0.454382},{-0.043758},{-0.840326},{-2.50552},

{0.681574},{0.615937},{-0.644552}}

− równania (10):

In[29]:= f0=Table[(x-x2[[1]])^k,{k,0,2}].{f[x2[[1]]],

D[f[x],x]/.x→x2[[1]], wsp[[1,1]]}//N

Out[29]:= 2.+0.454382x²

In[30]:= f1=Table[(x-x2[[2]])^k,{k,0,2}].{f[x2[[2]]],

wsp[[5,1]],wsp[[2,1]]}//N

Out[30]:= 2.25559+0.681574(-0.75+x)-0.043758(-0.75+x)²

In[31]:= f2=Table[(x-x2[[3]])^k,{k,0,2}].{f[x2[[3]]],

wsp[[6,1]],wsp[[3,1]]}//N

Out[31]:= 2.74216+0.615937(-1.5+x)-0.840326(-1.5+x)²

In[32]:= f3=Table[(x-x2[[4]])^k,{k,0,2}].{f[x2[[4]]],

wsp[[7,1]],wsp[[4,1]]}//N

Out[32]:= 2.73143-0.644552(-2.25+x)-2.50552(-2.25+x)²

− połączenie otrzymanych funkcji w całość:

In[33]:= int[x_]=If[x<xa[[2]],f0,If[x<xa[[3]],f1,If[x<xa[[4]],

f2,f3]]]

Out[33]:= If[x<0,f0,If[x<xa[[3]],f1,If[x<xa[[4]],f2,f3]]]

− wykres (linia gruba - interpolant, cienka - funkcja interpolowana):

In[34]:= Plot[{f[x],int[x]},{x,-1,3},GridLines→{x2,{0}},

PlotStyle→{Thickness[.003],Thickness[.006]}]

Out[34]:=

4. Rozwiązać równanie:
z warunkiem początkowym y(−1) = −2 w przedziale x ∈ (−1, 5) najpierw metodą Heuna z krokiem h = 0,2, potem analitycznie. Rozwiązania przedstawić na jednym wykresie.

Napisać równanie wielomianu Lagrange'a i Newtona, interpolującego funkcję y(x) w 4 węzłach rozłożonych wg Czebyszewa.

Po sprowadzeniu równania do postaci odpowiedniej dla obliczeń numerycznych mamy:
, a więc
. Zatem:

− pierwszy predyktor:

− pierwsze przybliżenie:

− korektor:
(12)

− drugi predyktor:

− właściwe przybliżenie - szukane rozwiązanie:

Wektor argumentów (przedział (−1, 5), skok h = 0,2):

In[1]:= x0 = Table[i,{i,−1,5,.2}]

Out[1]:= {-1.,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.,1.2,1.4,

1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.2,3.4,3.6,3.8,4.,4.2,

4.4,4.6,4.8,5.}

Szablony dla poszczególnych wielkości z algorytmu, czyli p1, y0, k, p2, y1:

In[2]:= y1=Table[0,{i,1,Length[x0]}]

Out[2]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0}

In[3]:= p1=Table[0,{i,1,Length[x0]}]

Out[3]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0}

In[4]:= y10=Table[0,{i,1,Length[x0]}]

Out[4]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0}

In[5]:= k=Table[0,{i,1,Length[x0]}]

Out[5]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0}

In[6]:= p2=Table[0,{i,1,Length[x0]}]

Out[6]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0}

Warunek początkowy:

In[7]:= y1[[1]] = -2

Out[7]:= -2

Pętla:

In[8]:= Do[p1[[i]]=0.5 Sin[2x0[[i]]]-y1[[i]]Cos[x0[[i]]];

y10[[i]]=y1[[i]]+.2 p1[[i]]; k[[i]]=0.5 Sin[2x0[[i+1]]]

-y10[[i]]Cos[x0[[i+1]]; p2[[i]]=
;

y1[[i+1]]=y1[[i]]+.2 p2[[i]],{i,1,Length[x0]-1}]

Wynik:

In[9]:= y1

Out[9]:= {-2,-1.85676,-1.68384,-1.48926,-1.28133,-1.06831,

-0.858247,-0.658826,-0.477236,-0.32,-0.192782,-0.100192,

-0.0456077,-0.0310326,-0.0570098,-0.122593,-0.225384,

-0.361629,-0.52637,-0.713629,-0.916615,-1.12791,

-1.33967,-1.54374,-1.73182,-1.89572,-2.02767,-2.12091,

-2.17029,-2.17298,-2.12883}

Rozwiązanie analityczne:

In[10]:= DSolve[{y'[x]+y[x]Cos[x]
Sin[2x],y[-1]-2},y[x],x]

Out[10]:= {{y[x]→e^{-Sin[1]-Sin[x]}(-1-e^{Sin[1]+Sin[x]}+Sin[1]+e^{Sin[1]+Sin[x]}

Sin[x])}}

In[11]:= f[x_]=Simplify[e^{-Sin[1]-Sin[x]}(-1-e^{Sin[1]+Sin[x]}+Sin[1]+

e^{Sin[1]+Sin[x]}Sin[x])]

Out[11]:= e^{-Sin[1]-Sin[x]}(-1-e^{Sin[1]+Sin[x]}+Sin[1])+Sin[x]

Wykresy:

In[12]:= Plot[f[x],{x,-1,5},Epilog→{Table[Point[{x0[[i]],

y1[[i]]}],{i,1,Length[x0]}]}]

Out[12]:=

Przybliżenie funkcji f(x) (wczytanej w instrukcji 11) wielomianem Lagrange'a w 4 węzłach (czyli wielomian 3. stopnia), rozłożonych wg Czebyszewa:

− węzły:
, gdzie i = 0, …, N - numer węzła, N - stopień wielomianu przybliżającego, a - początek przedziału, b - koniec przedziału (u nas N = 3,

a = −1, b = 5):

In[11]:= xc = Table
//N

Out[11]:= {-0.771639,0.85195,3.14805,4.77164}

− wielomian Lagrange'a: zgodnie z (6) i zasadą dotyczącą indeksowania (od 1, nie od 0)

⇒

In[12]:=

Out[12]:= 0.052004(-4.77164+x)(-3.14805+x)(-0.85195+x)-0.0191368

(-4.77164+x)(-3.14805+x)(0.771639 +x)+0.0735845

(-4.77164+x)(-0.85195+x)(0.771639+x)-0.0619007

(-3.14805+x)(-0.85195+x) (0.771639+x)

W postaci sumy:

In[13]:= Expand[WL]

Out[13]:= -0.784614+0.966407x-0.476561x²+0.0445511x³

Przybliżenie funkcji f(x) wielomianem Newtona w 4 węzłach - tych samych, jak w instrukcji 11:

− wielomian Newtona zgodnie z (7) i zasadą dotyczącą indeksowania

⇒
(13)

gdzie
, m = 0,…,3 ⇒

⇒
, m = 1,…,4

− wyrazy b: szablon będzie zawierał 4 wyrazy, bo we wz. (13) n zmienia się od 1 do N + 1 = 4

In[14]:= b = Table[0,{i,1,4}]

Out[14]:= {0,0,0,0,}

In[15]:= Do[
,{m,1,4}]

In[16]:= b

Out[16]:= {-1.83456,0.957709,-0.332734,0.0445511}

− wielomian:

In[17]:=

Out[17]:= -1.83456+0.957709(0.771639+x)-0.332734(-0.85195+x)

(0.771639+x)+0.0445511(-3.14805+x)(-0.85195+x)

(0.771639+x)

W postaci sumy - wzór jest ten sam, co dla WL:

In[18]:= Expand[WN]

Out[18]:= -0.784614+0.966407x-0.476561x²+0.0445511x³

Wykres (funkcja f(x) i wielomian Lagrange'a):

In[19]:= Plot[{f[x],WL},{x,-1,5},PlotStyle→{Thickness[.004],

Thickness[.008]}]

Out[19]:=

Przybliżenie jest dość słabe - trzeba zwiększyć liczbę węzłów (wystarczy o 1; jaki będzie wynik - sprawdźcie sami).

(x1[[2]]−x1[[1]])³

0

(x1[[3]]−x1[[2]])³

0

0

0

0

0

0

(x1[[4]]−x1[[3]])³

(x1[[5]]−x1[[4]])³

0

0

0

0

0

3(x1[[2]]−x1[[1]])²

0

0

0

0

0

0

0

0

0

3(x1[[3]]−x1[[2]])²

0

0

6(x1[[2]]−x1[[1]])

0

6(x1[[3]]−x1[[2]])

0

0

0

0

0

3(x1[[4]]−x1[[3]])²

6(x1[[4]]−x1[[3]])

0

0

0

0

3(x1[[5]]−x1[[4]])²

0

0

0

2(x1[[2]]−x1[[1]])

0

0

0

(x1[[2]]−x1[[1]])²

0

0

2

−2

0

0

2(x1[[3]]−x1[[2]])

0

0

0

(x1[[3]]−x1[[2]])²

0

0

2

-2

0

0

2(x1[[4]]−x1[[3]])

0

0

0

0

0

(x1[[4]]−x1[[3]])²

2

-2

0

0

2(x1[[5]]−x1[[4]])

0

0

0

0

0

0

(x1[[5]]−x1[[4]])²

0

0

0

0

0

-1

0

0

x1[[3]]−x1[[2]]

0

1

0

0

1

0

0

0

-1

0

0

x1[[4]]−x1[[3]]

0

0

0

0

1

0

0

-1

0

0

x1[[5]]−x1[[4]]

0

A=

(x2[[2]]−x2[[1]])³

0

(x2[[3]]−x2[[2]])³

0

0

0

(x2[[4]]−x2[[3]])³

(x2[[5]]−x2[[4]])³

0

0

0

0

0

0

0

0

3(x2[[2]]−x2[[1]])²

0

0

0

0

0

0

3(x2[[3]]−x2[[2]])²

0

0

0

6(x2[[2]]−x2[[1]])

0

0

6(x2[[3]]−x2[[2]])

0

0

0

3(x2[[4]]−x2[[3]])²

0

0

0

6(x2[[4]]−x2[[3]])

0

3(x2[[5]]−x2[[4]])²

0

0

0

(x2[[2]]−x2[[1]])²

0

0

0

2(x2[[2]]−x2[[1]])

0

0

0

2

0

0

0

(x2[[3]]−x2[[2]])²

0

0

0

2(x2[[3]]−x2[[2]])

0

0

−2

2

0

(x2[[4]]−x2[[3]])²

0

0

0

0

0

2(x2[[4]]−x2[[3]])

0

0

-2

2

(x2[[5]]−x2[[4]])²

0

0

0

0

0

0

2(x2[[5]]−x2[[4]])

0

0

-2

0

x2[[3]]−x2[[2]]

0

0

-1

0

0

0

1

0

0

0

x2[[4]]−x2[[3]]

0

0

-1

0

0

1

0

0

0

0

x2[[5]]−x2[[4]]

0

0

-1

0

0

1

0

0

0

A=

Wykreślanie wierszy i kolumn dla splajnu 3. stopnia w 4 węzłach

Wykreślanie wierszy i kolumn dla splajnu 2. stopnia w 5 węzłach

f[x1[[2]]]-f[x1[[1]]]-(D[f[x],x]/.x→x1[[1]])(x1[[2]]-x1[[1]])

f[x1[[3]]]-f[x1[[2]]]

f[x1[[4]]]-f[x1[[3]]]

-D[f[x],x]/.x→x1[[1]]

f[x1[[5]]]-f[x1[[4]]]

0

0

D[f[x],x]/.x→x1[[5]]

0

0

0

w=

4

w=

f[x2[[2]]]-f[x2[[1]]]-(D[f[x],x]/.x→x2[[1]])(x2[[2]]-x2[[1]])

f[x2[[3]]]-f[x2[[2]]]

f[x2[[4]]]-f[x2[[3]]]

f[x2[[5]]]-f[x2[[4]]]

-D[f[x],x]/.x→x2[[1]]

0

0

D[f[x],x]/.x→x2[[5]]

0

0

0

---