> restart:
> with(plots):
Interpolacja
Dane:
Zadanie funkcji dyskretnej
> 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=20,view=[-1..6,0..10],a xes=framed);
> pkt:=%:
Interpolacja za pomocą jednego wielomianu 1. Generowanie i rozwiązywanie układu równań
> w:=add(a[i]*x^i,i=0..n);
:=
2
3
4
5
w
a
a
a
0
a 1 x a 2 x
3 x
4 x
a 5 x
n
(
w x )
i
a x y , k 0.. n k
i
k
k
i0
> for i from 0 to n do
r||i:=eval(w,x=X[i])=Y[i];
end;
r0 := a 5
0
:=
r1
a 6
0
a 1 a 2 a 3 a 4 a 5
:=
r2
a 2 a 4 a 8 a 16 a 32 a 6
0
1
2
3
4
5
:=
r3
a 3 a 9 a 27 a 81 a 243 a 5
0
1
2
3
4
5
:=
r4
a 4 a 16 a 64 a 256 a 1024 a 4
0
1
2
3
4
5
:=
r5
a 5 a 25 a 125 a 625 a 3125 a 6
0
1
2
3
4
5
> roz:=solve({r||(0..n)},{seq(a[i],i=0..n)}); 29
-11
1
-1
1
:=
roz
{ a 5, a , a
, a , a , a
}
0
1
20 2
24
3
24 4 24 5 120
> assign(roz):
> w;
29
11
1
1
1
5
2
3
4
5
x
x
x
x
x
20
24
24
24
120
> wp:=plot(w,x=0..5):
> display({pkt,wp});
> eval(w,x=3.5);
4.371093750
2. Wielomian interpolacyjny Lagrange'a
wielomianypomocnicze
n
( x x )
L ( x)
j
i
x x
j 0, j i (
)
i
j
> 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 ) ( x2 ) ( x3 ) (
x
4 ) ( x5 )
L :=
0
120
x (
x
2 ) ( x3 ) ( x4 ) (
x
5 )
L :=
1
24
x (
x
1 ) ( x3 ) ( x4 ) (
x
5 )
L :=
2
12
x (
x
1 ) ( x2 ) ( x4 ) (
x
5 )
L :=
3
12
x (
x
1 ) ( x2 ) ( x3 ) (
x
5 )
L :=
4
24
x (
x
1 ) ( x2 ) ( x3 ) (
x
4 )
L :=
5
120
n
(
w x) y L ( x) i
i
i0
> wLagr:=expand(add(Y[i]*L[i],i=0..n)); # bez rozwiazywania ukladu równań 29
11
1
1
1
:=
2
3
4
5
wLagr
5
x
x
x
x
x
20
24
24
24
120
> w;
29
11
1
1
1
5
2
3
4
5
x
x
x
x
x
20
24
24
24
120
3. Komenda w Maple'u
> interp(X,Y,x);
29
11
1
1
1
5
2
3
4
5
x
x
x
x
x
20
24
24
24
120
> interp(X,Y,3.5);
4.371093750
> a:='a';
:=
a
a
Krzywe sklejane trzeciego stopnia
1. Generowanie i rozwiązywanie układu równań
> for i to n do s[i]:=a[i]+b[i]*x+c[i]*x^2+d[i]*x^3 end do;
3
s := a
1
1
b 1 x c 1 x
d 1 x
2
3
s := a
2
2
b 2 x c 2 x
d 2 x
2
3
s := a
3
3
b 3 x c 3 x
d 3 x
2
3
s := a
4
4
b 4 x c 4 x
d 4 x
2
3
s := a
5
5
b 5 x c 5 x
d 5 x
Warunki przechodzenia krzywych przez węzły interpolacji s ( x ) y , i 1.. n i
i
i
s ( x ) y
, i 1.. n
i
i 1
i 1
> 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 := a 6
1
b 1 c 1 d 1
r6 := a 5
1
r2 := a 2 b 4 c 8 d 6
2
2
2
2
r7 := a 6
2
b 2 c 2 d 2
:=
r3
a 3 b 9 c 27 d 5
3
3
3
3
r8 := a 2 b 4 c 8 d 6
3
3
3
3
r4 := a 4 b 16 c 64 d 4
4
4
4
4
:=
r9
a 3 b 9 c 27 d 5
4
4
4
4
:=
r5
a 5 b 25 c 125 d 6
5
5
5
5
:=
r10
a 4 b 16 c 64 d 4
5
5
5
5
Warunki równości pierwszych i drugich pochodnych w wewnętrznych węzłach interpolacji s( x ) s ( x ), i 1.. n 1
i
i
i 1
i
s (
x ) s ( x ), i 1.. n 1
i
i
i 1
i
> 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
b 2 c 3 d b 2 c 3
1
1
1
2
2
d 2
:=
r16
2 c 6 d 2 c 6
1
1
2
d 2
:=
r12
b 4 c 12 d b 4 c 12
2
2
2
3
3
d 3
:=
r17
2 c 12 d 2 c 12
2
2
3
d 3
:=
r13
b 6 c 27 d b 6 c 27
3
3
3
4
4
d 4
:=
r18
2 c 18 d 2 c 18
3
3
4
d 4
:=
r14
b 8 c 48 d b 8 c 48
4
4
4
5
5
d 5
:=
r19
2 c 24 d 2 c 24
4
4
5
d 5
Zerowanie drugich pochodnych w skrajnych węzłach s (
x ) s( x ) 0
1
0
n
n
> 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 c 0
1
r20 := 2 c 30 d 0
5
5
Rozwizanie układu równań
> 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)});
90
98
-388
1724
23
26
512
-1072
-15
-9
roz := { a 5, a
, a , a
, a
, b , b 2, b , b
, b
, c 0, c
, c ,
1
2
19 3 19 4
19
5
19
1
19 2
3
19 4
19
5
19
1
2
19
3
19
225
-4
1
18
-15
c -9, c
, d , d , d 0, d , d
}
4
5
19
1
19 2 19 3
4
19 5
19
> assign(roz):
> for i to n do s[i]; end do;
23
4
5
3
x
x
19
19
90
15
1
2
2
3
x
x
x
19
19
19
98 26
9
2
x
x
19 19
19
388 512
18
2
3
x9 x
x
19
19
19
1724 1072
225
15
2
3
x
x
x
19
19
19
19
Zapis funkcji przedziałami zmiennej
> S:=piecewise(x<X[1],s[1],x<X[2],s[2],x<X[3],s[3],x<X[4],s[4],s[5]);
23
4 3
5
x
x
x1
19
19
90
15
1
2
2
3
x
x
x2
x
19
19
19
98 26
9
:=
2
S
x
x
x3
19 19
19
388 512
18
2
3
x9 x
x
x
4
19
19
19
1724 1072 225
15
2
3
x
x
x
otherwise
19
19
19
19
> Sp:=plot(S,x=X[0]..X[n],color=blue):
> display({pkt,wp,Sp});
>
2. Komenda Maple'a (Spline w pakiecie CurveFitting)
> with(CurveFitting);
[ ArrayInterpolation, BSpline, BSplineCurve, Interactive, LeastSquares, PolynomialInterpolation, RationalInterpolation, Spline, ThieleInterpolation ]
> type(X,list);
false
> X:=convert(X,list): Y:=convert(Y,list):
> type(X,list);
true
> SM:=Spline(X,Y,x);
23
4 3
5
x
x
x1
19
19
90
15
1
2
2
3
x
x
x2
x
19
19
19
98 26
9 2
SM :=
x
x
x3
19 19
19
388 512
18
2
3
x
9 x
x
x
4
19
19
19
1724 1072 225
15
2
3
x
x
x
otherwise
19
19
19
19
> SM,S;
23
4
23
4
3
3
5
x
x
x1
5
x
x
x
1
19
19
19
19
90
15
1
90
15
1
2
2
3
2
3
x
x
x2
2 x
x
x
x
x
2
19
19
19
19
19
19
98 26
9
98 26
9
2
2
x
x
x3
,
x
x
x
3
19 19
19
19 19
19
388 512
18
388 512
18
2
3
2
3
x
9 x
x
x9 x
x
4
x
x
4
19
19
19
19
19
19
1724 1072
225
15
1724 1072
225
15
2
3
2
3
x
x
x
x
x
otherwise
x
otherwise
19
19
19
19
19
19
19
19
>