C4MW.html

多項式補間(polynomial interpolation)

> restart;
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]);

> with(plots):

Warning, the name changecoords has been redefined

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

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

> with(LinearAlgebra):

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

> display(f1p,l1p);

> 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;

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

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

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

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

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

f20の取るべき値の導出

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

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

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

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

> evalb(fac_f2p=fac_newton);

Numerical integration

> restart;
f1:=x->4/(1+x^2);
plot(f1(x),x=0..5);

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

> int(1/(1+x^2),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);

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);

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);

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)));

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

> evalb(e1=e2);