C3MW.html

FindRootAll.mw

fsolve

次のような二次方程式
x*x-4.0*x+1=0
の根を求める．先ず，関数の概形を表示してみる．

> plot(x*x-4.0*x+1,x=-1..5);

Mapleで根を求めるには，solveを用いる．

> (xp,xm):=solve(x^2-4*x+1=0);

> evalf(xm);

> fsolve(x*x-4.0*x+1=0);

ニュートン法の視覚化

> with(Student[Calculus1]):
NewtonsMethodTutor(x^2-4*x+1, x=-1..2);

Schematic illustrations for bisection method

> restart:
with(plots):
with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

> f1:=x->x^2-4*x+1;

> p1:=plot(f1(x),x=0..2):

> c0:=pointplot([0,1],symbolsize=20,symbol=circle):
c1:=pointplot([2,f1(2)],symbolsize=20,symbol=circle):
c2:=pointplot([1,0],symbolsize=20,symbol=circle):
c3:=pointplot([1,f1(1)],symbolsize=20,symbol=circle):
c4:=pointplot([0.5,0],symbolsize=20,symbol=circle):
c5:=pointplot([0.5,f1(0.5)],symbolsize=20,symbol=circle):
l0:=arrow([0,f1(0)],[0,0],.05,.1,.2,color=green):
l1:=arrow([2,f1(2)],[2,0],.05,.1,.2,color=green):
l2:=arrow([1,0],[1,f1(1)],.05,.1,.2,color=green):
l3:=arrow([0.5,0],[0.5,f1(0.5)],.05,.1,.2,color=green):

> d1:=display(p1):
d2:=display(p1,c0,c1):
d3:=display(l0,l1,p1,c0,c1,c2):
d4:=display(l0,l1,p1,c0,c1,c2,c3,l2):
d5:=display(l0,p1,c0,c2,c3,l2,c4,c5,l3):

> display([d1,d2,d3,d4,d5],insequence=true,view=[-0.2..2.2,-3.2..1.2]);

> display([d1,d2,d3,d4,d5],view=[-0.2..2.2,-3.2..1.2]);

Schematic illustrations for Newton-Raphson method

> restart:
with(plots):
with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

> f1:=x->x^2-4*x+1;
df1:=(x,x0)->(2*x0-4)*(x-x0)+f1(x0);

> df1(x,0);

> p1:=plot({f1(x)},x=-0.5..2):
p2:=plot({df1(x,1),f1(x)},x=-0.5..2):
p3:=plot({df1(x,0),f1(x)},x=-0.5..2):

> c0:=pointplot([1,f1(1)],symbolsize=20,symbol=circle):
a1:=arrow([0,0],[0,f1(0)],.05,.1,.2,color=green):
c1:=pointplot([0,f1(0)],symbolsize=20,symbol=circle):
a2:=arrow([0.25,0],[0.25,f1(0.25)],.02,.04,.01,color=green):
c2:=pointplot([0.25,f1(0.25)],symbolsize=20,symbol=circle):
d1:=display(p1):
d2:=display(p1,c0):
d3:=display(p2,c0):
d35:=display(p2,c0,a1):
d4:=display(p2,c1):
d5:=display(p3,a2):
d6:=display(p3,a2,c2):

> display([d1,d2,d3,d35,d4,d5,d6],insequence=true,view=[-0.1..1.1,-2.1..1.1]);

> display([d1,d2,d3,d35,d4,d5,d6],view=[-0.1..1.1,-2.1..1.1]);

Derivation of eq.(3.11) on p.19

解の収束の速さ

> restart;
#for bisection method
e[n+1]=e[n]/2;

> s1:=series(f(xp-e),e,2);

> ee:=solve(convert(s1,polynom)=0,e);

> xpp:=xp-ee;

> f:=x->a*(x-x0)*(x-x1);

> xpp:=xp-f(xp)/subs(x=xp,diff(f(x),x));

> xpp2:=subs({xp-x1=dx1,xp-x0=dx0},xpp);

> simplify(subs(dx1=dx0/alpha,xpp2-xp));

> series(%,alpha,3);

> subs({dx0=xp-x0,dx1=xp-x1},subs(alpha=dx0/dx1,series(simplify(subs(dx1=dx0/alpha,xpp2)),alpha,3)));

Bairstow-Hitchcock method

Simplest explanation

> restart;

> P:=sort(sum('x^i*a[i]','i'=0..n),x);

> n:=3;
a:=array(0..n+1,[-40,62,-23,1]):
sort(P,x);

> factor(P);

> plot(P,x=-10..30);
plot(P,x=-1..21,y=-10..10);

> fsolve(P=0,x);

Analytical solution

Bairstow-Hitchcock法の基本的な考え方の導出．

> restart:
n:=3;

> P1:=sort(sum(x^i*a[i],i=0..n),x);

> Q:=sum(x^i*b[i],i=0..n-2);

> P2:=sort(expand((x^2+B*x+C)*Q+R*x+S),x);

二つの式f1とf2を等値して，係数が満たすべき方程式を導出．

> for i from 0 to n do
eq||i:=coeff(P1,x,i)=coeff(P2,x,i);
end do;

係数の下の方のbの満たすべき式を導出．

> b[n-2]:=1;
for i from n-3 to 0 by -1 do
b[i]:=solve(eq||(i+2),b[i]);
end do;

r,sが０となる方程式をたてる．

> Eq_R:=solve(expand(eq1),R)=0;
Eq_S:=solve(expand(eq0),S)=0;

> a:=array(0..n,[-40,62,-23,1]);

上のような係数を与えて，p,qを求める．

> fsolve({Eq_R,Eq_S},{B,C});

Numerical solution

数値解を求めるコード

初期値の代入．係数の添え字が逆になっている事に注意．

> restart;
n:=3;
a:=array(0..n+2,[0,1,-23,62,-40]):
b:=array(0..n+2,[0,1,0,0,0,0]):
d:=array(0..n+2,[0,1,0,0,0,0]):

> B:=1.0;
C:=1.0;

以下の二つのコマンド領域を収束するまで繰り返す．

> for i from 2 to n+2 do
j:=i-1;
k:=i-2;
b[i]:=a[i]-B*b[j]-C*b[k];
d[i]:=b[i]-B*d[j]-C*d[k];
od:

> j:=n-1;
k:=n-2;
det1:=d[j]*d[j]-(d[n]-b[n])*d[k];
dB:=(b[n]*d[j]-b[n+1]*d[k])/det1;
dC:=(b[n+1]*d[j]-(d[n]-b[n])*b[n])/det1;
B:=B+dB;
C:=C+dC;

> print(convert(b,list));

Note for Numerical solution

det1, dB, とdCはテキスト(3.41)式の逆行列から出てきた式．

> with(LinearAlgebra):

> A:=Matrix(1..2,1..2,[[a11,a12],[a21,a22]]);

> MatrixInverse(A);

ちょっと込み入った例をMapleで扱うと．．．

> restart:
n:=4;
P:=sort(sum('x^i*a[i]','i'=0..n),x);
a:=array(0..n+1,[-41,62,-23,1,1]):
sort(f1,x);

> fsolve(P=0,x,complex);

> x1,x2:=fsolve(P=0,x);

> Q:=x^2-(x1+x2)*x+x1*x2;

> R:=factor(P/Q);

> expand(R*Q);