> restart:
> with(plots):with(LinearAlgebra):
Interpolacja
Dane:
> ?VandermondeMatrix
> n:=5;
n := 5
> X:=array(0..n,[0,1,2,3,4,5]):
> Y:=array(0..n,[5,6,6,5,4,6]):
> plot(zip((x,y)->[x,y],X,Y),style=point,symbol=cross,symbolsize=2
0,view=[-1..6,0..10],axes=framed);
> pkt:=%:
Jawne generowanie i rozwiazywanie układu równan
> w:=add(a[i]*x^i,i=0..n);
w := a0 + a1 x + a2 x2 + a3 x3 + a4 x4 + a5 x5
> for i from 0 to n do
r||i:=eval(w,x=X[i])=Y[i];
end;
r0 := a0 = 5
r1 := a0 + a1 + a2 + a3 + a4 + a5 = 6
r2 := a0 + 2 a1 + 4 a2 + 8 a3 + 16 a4 + 32 a5 = 6
r3 := a0 + 3 a1 + 9 a2 + 27 a3 + 81 a4 + 243 a5 = 5
r4 := a0 + 4 a1 + 16 a2 + 64 a3 + 256 a4 + 1024 a5 = 4
r5 := a0 + 5 a1 + 25 a2 + 125 a3 + 625 a4 + 3125 a5 = 6
> roz:=solve({r||(0..n)},{seq(a[i],i=0..n)});
1 -1 1 -11 29
roz := { a0 = 5, a5 = , a4 = , a3 = , a2 = , a1 = }
120 24 24 24 20
> assign(roz):
> w;
29 11 1 1 1
5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
> wr:=copy(w); # współczynniki a[i] ulegną zmianie podczas
dalszych obliczeń
29 11 1 1 1
wr := 5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
> pwr:=plot(wr,x=X[0]..X[n]):
> display([pkt,pwr]);
Macierz Van der Monde'a
Ai+1, j+1 = xi j , i, j = 0..n
bi+1 = yi , i = 0..n
> A:=Matrix(n+1):b:=Vector(n+1):
> for i from 0 to n do
b[i+1]:=Y[i]:
for j from 0 to n do
A[i+1,j+1]:=X[i]^j:
end do:
end do:
> A,b;
îÅ‚1 0 0 0 0 0Å‚Å‚ îÅ‚ 5Å‚Å‚
ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł
6śł
ïÅ‚1 1 1 1 1 1śł ïÅ‚
ïÅ‚ śł ïÅ‚ śł
6śł
ïÅ‚1 2 4 8 16 32śł, ïÅ‚
ïÅ‚1 3 9 27 81 243śł ïÅ‚
5śł
ïÅ‚ śł ïÅ‚ śł
ïÅ‚1 4 16 64 256 1024śł ïÅ‚
4śł
ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł
ðÅ‚1 5 25 125 625 3125ûÅ‚ ðÅ‚ 6ûÅ‚
> a:=A^(-1).b;
îÅ‚ 5 Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
29
ïÅ‚ śł
ïÅ‚ śł
20
ïÅ‚ śł
ïÅ‚ śł
-11
ïÅ‚ śł
ïÅ‚ śł
ïÅ‚ śł
24
ïÅ‚ śł
a := ïÅ‚ 1 śł
ïÅ‚ śł
ïÅ‚ śł
24
ïÅ‚ śł
ïÅ‚ śł
-1
ïÅ‚ śł
ïÅ‚ śł
ïÅ‚ śł
24
ïÅ‚ śł
ïÅ‚ śł
1
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚ 120 ûÅ‚
> wV:=add(a[i+1]*x^i,i=0..n); # wielomian interpolacyjny
wyznaczony z ukadu równan A.a=b
29 11 1 1 1
wV := 5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
Wielomian interpolacyjny Lagrange'a
wielomianypomocnicze
n
(x - xj )
Li (x) =
"
(xi - xj )
j=0, j`"i
> for i from 0 to n do
L[i]:=factor(mul((x-X[j])/(X[i]-X[j]),j=0..i-1)*mul((x-X[j])/(X[
i]-X[j]),j=i+1..n)):
end do;
( x - 1) ( x - 2 ) ( x - 3) ( x - 4 ) ( x - 5 )
L0 := -
120
x ( x - 2) (x - 3 ) ( x - 4 ) (x - 5 )
L1 :=
24
x ( x - 1) (x - 3 ) ( x - 4 ) (x - 5 )
L2 := -
12
x ( x - 1) (x - 2 ) ( x - 4 ) (x - 5 )
L3 :=
12
x ( x - 1) (x - 2 ) ( x - 3 ) (x - 5 )
L4 := -
24
x ( x - 1) (x - 2 ) ( x - 3 ) (x - 4 )
L5 :=
120
n
w(x) = yi Li (x)
"
i=0
> wL:=expand(add(Y[i]*L[i],i=0..n)); # bez rozwiazywania ukladu
równań
29 11 1 1 1
wL := 5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
Komenda Maple'a (interp)
> wM:=interp(X,Y,x);
29 11 1 1 1
wM := 5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
Porównanie
> wr;wV;wL;wM;
29 11 1 1 1
5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
29 11 1 1 1
5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
29 11 1 1 1
5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
29 11 1 1 1
5 + x - x2 + x3 - x4 + x5
20 24 24 24 120
> a:='a':b:='b':
Krzywe sklejane trzeciego stopnia - generowanie i rozwiazywanie ukladu równan
> for i to n do s[i]:=a[i]+b[i]*x+c[i]*x^2+d[i]*x^3 end do;
s1 := a1 + b1 x + c1 x2 + d1 x3
s2 := a2 + b2 x + c2 x2 + d2 x3
s3 := a3 + b3 x + c3 x2 + d3 x3
s4 := a4 + b4 x + c4 x2 + d4 x3
s5 := a5 + b5 x + c5 x2 + d5 x3
Warunki przechodzenia krzywych przez wezły interpolacji
> for i to n do
r||i:=eval(s[i],x=X[i])=Y[i];
r||(n+i):=eval(s[i],x=X[i-1])=Y[i-1];
end do;
r1 := a1 + b1 + c1 + d1 = 6
r6 := a1 = 5
r2 := a2 + 2 b2 + 4 c2 + 8 d2 = 6
r7 := a2 + b2 + c2 + d2 = 6
r3 := a3 + 3 b3 + 9 c3 + 27 d3 = 5
r8 := a3 + 2 b3 + 4 c3 + 8 d3 = 6
r4 := a4 + 4 b4 + 16 c4 + 64 d4 = 4
r9 := a4 + 3 b4 + 9 c4 + 27 d4 = 5
r5 := a5 + 5 b5 + 25 c5 + 125 d5 = 6
r10 := a5 + 4 b5 + 16 c5 + 64 d5 = 4
Warunki równości pierwszych i drugich pochodnych w wewnetrznych wezlach interpolacji
> for i to n-1 do
r||(2*n+i):=eval(diff(s[i],x),x=X[i])=eval(diff(s[(i+1)],x),x=X[
i]);
r||(3*n+i):=eval(diff(s[i],x,x),x=X[i])=eval(diff(s[(i+1)],x,x),
x=X[i]);
end do;
r11 := b1 + 2 c1 + 3 d1 = b2 + 2 c2 + 3 d2
r16 := 2 c1 + 6 d1 = 2 c2 + 6 d2
r12 := b2 + 4 c2 + 12 d2 = b3 + 4 c3 + 12 d3
r17 := 2 c2 + 12 d2 = 2 c3 + 12 d3
r13 := b3 + 6 c3 + 27 d3 = b4 + 6 c4 + 27 d4
r18 := 2 c3 + 18 d3 = 2 c4 + 18 d4
r14 := b4 + 8 c4 + 48 d4 = b5 + 8 c5 + 48 d5
r19 := 2 c4 + 24 d4 = 2 c5 + 24 d5
Zerowanie drugich pochodnych w skrajnych węzlach
> r||(3*n):=eval(diff(s[1],x,x),x=X[0])=0;
r||(4*n):=eval(diff(s[n],x,x),x=X[n])=0;
r15 := 2 c1 = 0
r20 := 2 c5 + 30 d5 = 0
> roz:=solve({r||(1..4*n)},{seq(a[i],i=1..n),seq(b[i],i=1..n),seq(
c[i],i=1..n),seq(d[i],i=1..n)});
-1072 -15 225 23 98 -4 26
roz := { b5 = , d5 = , c5 = , b1 = , a1 = 5, c1 = 0, a3 = , d1 = , b3 = , d3 = 0,
19 19 19 19 19 19 19
18 -388 512 1724 -9 1 90 -15
d4 = , a4 = , b4 = , a5 = , c3 = , d2 = , a2 = , b2 = 2, c4 = -9, c2 = }
19 19 19 19 19 19 19 19
> assign(roz):
> for i to n do s[i]; end do;
23 4
5 + x - x3
19 19
90 15 1
+ 2 x - x2 + x3
19 19 19
98 26 9
+ x - x2
19 19 19
388 512 18
- + x - 9 x2 + x3
19 19 19
1724 1072 225 15
- x + x2 - x3
19 19 19 19
Wygenerowanie jednej funkcji skladajacej sie z poszczególnych funkcji sklejanych
> S:=piecewise(x]);
Å„Å‚ 23 4
ôÅ‚
5 + x - x3 x < 1
ôÅ‚
19 19
ôÅ‚
ôÅ‚
90 15 1
ôÅ‚
+ 2 x - x2 + x3 x < 2
ôÅ‚
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚ 98 26 9
S := ôÅ‚ + x - x2 x < 3
òÅ‚
ôÅ‚
19 19 19
ôÅ‚
ôÅ‚
388 512 18
ôÅ‚
ôÅ‚ - + x - 9 x2 + x3 x < 4
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚
1724 1072 225 15
ôÅ‚
- x + x2 - x3 otherwise
ôÅ‚
ół 19 19 19 19
> pS:=plot(S,x=X[0]..X[n],color=blue):
> display({pkt,pwr,pS});
>
Komenda Maple'a (Spline w pakiecie CurveFitting)
> with(CurveFitting);
Warning, the name LeastSquares has been rebound
[BSpline, BSplineCurve, Interactive, LeastSquares, PolynomialInterpolation,
RationalInterpolation, Spline, ThieleInterpolation ]
> type(X,list);
false
> X:=convert(X,list);Y:=convert(Y,list);
X := [0, 1, 2, 3, 4, 5 ]
Y := [5, 6, 6, 5, 4, 6 ]
> type(X,list);
true
> SM:=Spline(X,Y,x,degree=3);
Å„Å‚ 23 4
ôÅ‚
5 + x - x3 x < 1
ôÅ‚
19 19
ôÅ‚
ôÅ‚
90 15 1
ôÅ‚
+ 2 x - x2 + x3 x < 2
ôÅ‚
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚ 98 26 9
SM := ôÅ‚ + x - x2 x < 3
òÅ‚
ôÅ‚
19 19 19
ôÅ‚
ôÅ‚
388 512 18
ôÅ‚
ôÅ‚ - + x - 9 x2 + x3 x < 4
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚
1724 1072 225 15
ôÅ‚
- x + x2 - x3 otherwise
ôÅ‚
ół 19 19 19 19
> S,SM;
Å„Å‚ 23 4
ôÅ‚
5 + x - x3 x < 1
ôÅ‚
19 19
ôÅ‚
ôÅ‚
90 15 1
ôÅ‚
+ 2 x - x2 + x3 x < 2
ôÅ‚
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚ 98 26 9
ôÅ‚
òÅ‚ + x - x2 x < 3 ,
ôÅ‚
19 19 19
ôÅ‚
ôÅ‚
388 512 18
ôÅ‚
ôÅ‚ - + x - 9 x2 + x3 x < 4
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚
1724 1072 225 15
ôÅ‚
- x + x2 - x3 otherwise
ôÅ‚
ół 19 19 19 19
Å„Å‚ 23 4
ôÅ‚
5 + x - x3 x < 1
ôÅ‚
19 19
ôÅ‚
ôÅ‚
90 15 1
ôÅ‚
+ 2 x - x2 + x3 x < 2
ôÅ‚
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚ 98 26 9
ôÅ‚
òÅ‚ + x - x2 x < 3
ôÅ‚
19 19 19
ôÅ‚
ôÅ‚
388 512 18
ôÅ‚
ôÅ‚ - + x - 9 x2 + x3 x < 4
ôÅ‚ 19 19 19
ôÅ‚
ôÅ‚
1724 1072 225 15
ôÅ‚
- x + x2 - x3 otherwise
ôÅ‚
ół 19 19 19 19
Interpolacja wielomianem Taylora - interpolacja (linearyzacja) energii potencjalnej w polu
grawitacyjnym
>
(i)
k
f (x0)
w(x) = (x - x0)i
"
i!
i=0
> y:=x->exp(-x^2);
( -x2 )
y := x e
> x0:=0;
x0 := 0
> k:=10; # k:=50;
k := 10
> wT:=add((D@@i)(y)(x0)/i!*(x-x0)^i,i=0..k);
1 1 1 1
wT := 1 - x2 + x4 - x6 + x8 - x10
2 6 24 120
> wt:=convert(taylor(y(x),x,k+1),polynom); # inny sposób
wyznaczania szeregu Taylora (szybszy!)
1 1 1 1
wt := 1 - x2 + x4 - x6 + x8 - x10
2 6 24 120
Porównanie dokladności wielomianu Taylora
> evalf(y(2),20);evalf(eval(wt,x=2),20);
0.018315638888734180294
-3.5333333333333333333
> plot([y(x),wt],x=-2..2,-1..1);
Zwiekszyć stopnień wielomianu k = 50
>
>