Integral.mw

多項式補間(polynomial interpolation)

> restart;
X:=[0,1,2,3];

Y:=[1,2,3,-2];

X := [0, 1, 2, 3]

Y := [1, 2, 3, -2]

> with(stats):

> eq_fit:=fit[leastsquare[[x,y],y=a3*x^3+a2*x^2+a1*x+a0,
{a3,a2,a1,a0}]]([X,Y]);

eq_fit := y = -x^3+3*x^2-x+1

> with(plots):

Warning, the name changecoords has been redefined

> f1:=unapply(rhs(eq_fit),x);

f1 := proc (x) options operator, arrow; -x^3+3*x^2-x+1 end proc

> f1p:=plot(f1(x),x=0..3):

> with(LinearAlgebra):

> list1:=[X,Y];
l1p:=listplot(Transpose(Matrix(list1))):

list1 := [[0, 1, 2, 3], [1, 2, 3, -2]]

> display(f1p,l1p);

[Plot]

> A:=Matrix(4,4):
for i from 1 to 4 do

for j from 1 to 4 do

 A[i,j]:=X[i]^(j-1);

end do;

end do:

> A;

Matrix([[1, 0, 0, 0], [1, 1, 1, 1], [1, 2, 4, 8], [1, 3, 9, 27]])

> MatrixInverse(A).Transpose(Matrix(Y));

Matrix([[1], [-1], [3], [-1]])

ニュートンの補間公式の導出.

関数F(x)をxの多項式として展開.その時の,係数の取るべき
値と,差商公式で得られる値が一致.

> restart:
F:=x->f0+(x-x0)*f1p+(x-x0)*(x-x1)*f2p;

F := proc (x) options operator, arrow; f0+(x-x0)*f1p+(x-x0)*(x-x1)*f2p end proc

> F(x1);
s1:=solve(F(x1)=f1,f1p);

f0+(x1-x0)*f1p

s1 := (f0-f1)/(-x1+x0)

f20の取るべき値の導出

> s2:=solve(F(x2)=f2,f2p);
fac_f2p:=factor(subs(f1p=s1,s2));

s2 := -(f0+f1p*x2-f1p*x0-f2)/((x2-x0)*(x2-x1))

fac_f2p := -(-f0*x1+x2*f0-x2*f1+x0*f1+f2*x1-f2*x0)/((-x1+x0)*(x2-x0)*(x2-x1))

ニュートンの差分商公式を変形

> ff11:=(f0-f1)/(x0-x1);
ff12:=(f1-f2)/(x1-x2);

ff2:=(ff11-ff12)/(x0-x2);

fac_newton:=factor(ff2);

ff11 := (f0-f1)/(-x1+x0)

ff12 := (f1-f2)/(x1-x2)

ff2 := ((f0-f1)/(-x1+x0)-(f1-f2)/(x1-x2))/(-x2+x0)

fac_newton := -(-f0*x1+x2*f0-x2*f1+x0*f1+f2*x1-f2*x0)/((-x1+x0)*(x2-x0)*(x2-x1))

二式が等しいかどうかをevalbで判定

> evalb(fac_f2p=fac_newton);

true

Numerical integration

> restart;
f1:=x->4/(1+x^2);

plot(f1(x),x=0..5);

f1 := proc (x) options operator, arrow; 4/(1+x^2) end proc

[Plot]

> int(f1(x),x=0..1);

Pi

> int(1/(1+x^2),x);

arctan(x)

> N:=8:
x0:=0:

xn:=1:

Digits:=20:

Midpoint rule(中点法)

> h:=(xn-x0)/N:
S:=0:

for i from 0 to N-1 do

 xi:=x0+(i+1/2)*h;

 dS:=h*f1(xi);

 S:=S+dS;

end do:

evalf(S);

3.1428947295916887799

Trapezoidal rule(台形公式)

> h:=(xn-x0)/N:
S:=f1(x0)/2:

for i from 1 to N-1 do

 xi:=x0+i*h;

 dS:=f1(xi);

 S:=S+dS;

end do:

S:=S+f1(xn)/2:

evalf(h*S);

3.1389884944910890093

Simpson's rule(シンプソンの公式)

> M:=N/2:
h:=(xn-x0)/(2*M):

Seven:=0:

Sodd:=0:

for i from 1 to 2*M-1 by 2 do

 xi:=x0+i*h;

 Sodd:=Sodd+f1(xi);

end do:

for i from 2 to 2*M-1 by 2 do

 xi:=x0+i*h;

 Seven:=Seven+f1(xi);

end do:

evalf(h*(f1(x0)+4*Sodd+2*Seven+f1(xn))/3);

3.1415925024587069144

Derivation of Simpson's rule

> restart;
f:=x->a*x^2+b*x+c;

e1:=expand(subs(x1=x0+h,int(f(x),x=x0..x1)));

f := proc (x) options operator, arrow; a*x^2+b*x+c end proc

e1 := a*x0^2*h+a*x0*h^2+1/3*a*h^3+b*x0*h+1/2*b*h^2+c*h

> e2:=expand(h/6*(f(x0)+4*f(x0+h/2)+f(x0+h)));

e2 := a*x0^2*h+a*x0*h^2+1/3*a*h^3+b*x0*h+1/2*b*h^2+c*h

> evalb(e1=e2);

true