> restart:

> with(plots):

> with(CurveFitting);

[ ArrayInterpolation, BSpline, BSplineCurve, Interactive, LeastSquares, PolynomialInterpolation, RationalInterpolation, Spline, ThieleInterpolation ]

Aproksymacja

Aproksymacja funkcji dyskretnej wielomianami

> n:=10;

:=

n

10

> lx:=array(0..n,[0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]):

> ly:=array(0..n,[1.,0.97,0.84,0.7,0.51,0.38,0.22,0.15,0.07,0.04,0.02]):

> plot(zip((x,y)->[x,y],lx,ly),style=point,symbol=cross,symbolsize=20);

> pkt:=%:

> m:=4; Digits:=20; # dla m=10

m := 4

:=

Digits

20

Generowanie i rozwiązywanie układu równań na współczynniki wielomianu aproksymującego a[ j]

n

m

j 2

S  [ y   a x ]

k

j

k

k 0

j 0

> S:=add((ly[k]-add(a[j]*lx[k]^j,j=0..m))^2,k=0..n); 2

2

S := ( 1. a ) ( 0.971.0 a 0.2 a 0.04 a 0.008 a 0.0016 a ) 0

0

1

2

3

4

2

2

( 0.841.0 a 0.4 a 0.16 a 0.064 a 0.0256 a ) ( 0.71.0 a 0.6 a 0.36 a 0.216 a 0.1296 a ) 0

1

2

3

4

0

1

2

3

4

2

2

( 0.511.0 a 0.8 a 0.64 a 0.512 a 0.4096 a ) ( 0.381.0 a 1.0 a 1.00 a 1.000 a 1.0000 a ) 0

1

2

3

4

0

1

2

3

4

2

2

( 0.221.0 a 1.2 a 1.44 a 1.728 a 2.0736 a ) ( 0.151.0 a 1.4 a 1.96 a 2.744 a 3.8416 a ) 0

1

2

3

4

0

1

2

3

4

2

2

( 0.071.0 a 1.6 a 2.56 a 4.096 a 6.5536 a ) ( 0.041.0 a 1.8 a 3.24 a 5.832 a 10.4976 a ) 0

1

2

3

4

0

1

2

3

4

2

( 0.021.0 a 2.0 a 4.00 a 8.000 a 16.0000 a ) 0

1

2

3

4

> for i from 0 to m do

r||i:=diff(S,a[i])=0;

end do;

:=

r0

22.00 a 22.00 a 30.800 a 48.4000 a 81.06560 a 9.8000

0

1

2

3

4

:=

r1

22.00 a 30.80 a 48.400 a 81.0656 a 141.32800 a 4.8720

0

1

2

3

4

:=

r2

30.800 a 48.400 a 81.0656 a 141.32800 a 253.235840 a 4.26240

0

1

2

3

4

:=

r3

48.4000 a 81.0656 a 141.32800 a 253.235840 a 462.8588800 a 4.651200

0

1

2

3

4

:=

r4

81.06560 a 141.32800 a 253.235840 a 462.8588800 a 858.78442496 a 5.8675200

0

1

2

3

4

> roz:=solve({r||(0..m)},{seq(a[i],i=0..m)}); roz := { a 1.0017482517482517483, a 0.078331390831390831391, a -1.5270250582750582751, 0

1

2

a 1.0035450660450660451, a -0.19121503496503496503 }

3

4

> assign(roz):

> g:=add(a[j]*x^j,j=0..m); # wielomian aproksymujący 2

3

g := 1.00174825174825174830.078331390831390831391 x1.5270250582750582751 x 1.0035450660450660451 x 0.19121503496503496503 4



x

> p:=plot(g,x=0..2):

> display([pkt,p]);

> S;

0.00094871794871794871841

> delta:=sqrt(S);

:=



0.030801265375272307446

> a:='a':

Komenda w Maple'u

Wrócić do m=4

> X:=convert(lx,list): Y:=convert(ly,list):

> w:=add(a[i]*x^i,i=0..m);

:=

2

3

4

w

a 



 a



0

a 1 x a 2 x 3 x

a 4 x

> gM:=LeastSquares(X,Y,x,curve=w); 2

3

gM := 1.00174825174825174830.078331390831390832734 x1.5270250582750582784 x 1.0035450660450660475 x 0.19121503496503496561 4



x

> g;gM;

1.0017482517482517483 0.078331390831390831391

2

3



x1.5270250582750582751 x 1.0035450660450660451 x 0.19121503496503496503 4



x

1.0017482517482517483 0.078331390831390832734

2

3



x1.5270250582750582784 x 1.0035450660450660475 x 0.19121503496503496561 4



x

> a:='a': # zmienić wartość stopnia wielomianu "m"

>

> restart:

Aproksymacja funkcji okresowej wielomianami trygonometrycznymi

> f:=piecewise(x<-3,x+4,x<-1,-2-x,x<1,x,x<3,2-x,-4+x); # kilka okresow na wykresie

 x



4

x-3

 2 x

x-1

:=

f

 x

x1

 

2

3



x

x

 4 x

otherwise

> plot(f,x=-5..5,scaling=constrained);

>

T 2

2

 2 



( )cos

i

a

f x

x dx , i  0,1,...



i





T

 T



T

 2

T 2

2

 2 



( )sin

i

b

f x

x dx , i  1, 2,...



i





T

 T



T

 2

> T:=4;

:=

T

4

> m:=5; # 5, 10, 15, 25, 50

m := 5

> for i from 0 to m do

a[i]:=2/T*evalf(Int(f*cos(2*Pi/T*i*x),x=-T/2..T/2)): b[i]:=2/T*evalf(Int(f*sin(2*Pi/T*i*x),x=-T/2..T/2)): end do:

m

1

2 i



2

g( x)  a  ( a cos x  b sin i

 x)

2

0

i

T

i

T

i 1



> g:=1/2*a[0]+add(a[i]*cos(2*Pi/T*i*x)+b[i]*sin(2*Pi/T*i*x),i=1..m);

  x 

 3  x 

 5  x 

:=

g

0.8105694690



sin

0.09006327436 

sin

0.03242277876 

sin



 2 

 2 

 2 

> plot([f,g],x=-5..5,color=[red,blue],scaling=constrained,numpoints=200); # porównanie

>