Tutorial 3: introduction to programming and Newton Method

Loop （繰り返し）　and If文

Do[  Task(counter) , { counter , start , end   } ]

Means repeat the task depending on the integer counter from the value counter=start  to the value counter=end.
ということは、整数 counter に依存するTask(counter) をcounter=start　から　counter=end　まで繰り返すことである。

Remark:  Do[  Task(counter) , { counter , end   } ]   (without a value “start”) means from counter=1 to counter=end.

Example:

f[x_]:=x^2
Do[   Print["Hello"] ;  Print[f[i]]  ,  {i , 10}  ]

P:={} ;   Do[  P = P~Union~{i}   ,   {i,4,10}   ]

If[ condition  ,   then  ,  else ]        If [ 条件  , 　満たすなら　, 満たさないなら　]

Examples

If[1 == 2, True, False]    (*  Comment: If 1 equal 2 then return True else Return False *)

n=1/2
If[n > 0 && n < 1, Print["n is between 0 and 1"], Print["n is negative or greater than 1"]]

Boolean test: AND , OR

Example:  Is 1>0  AND  1<2 ?       Is     ( 2 equal 1*1 )  or  ( 2 equal 1+1 )  ?

1 >0 && 1<2
2== 1*1 || 2== 1+1

If[ condition 1 || condition 2 ,   then  ,  else ]        If [ 条件１  ||　　条件２   , 　満たすなら　, 満たさないなら　]　
If[ condition 1 && condition 2  ,   then  ,  else ]        If [ 条件１　&&　　条件２  , 　満たすなら　, 満たさないなら　]

Programming  Input （入力）　and Output (出力）

Example:

factorial[n_]:=(
m:=1;
Do[m=m*i,{i,n}];
Return[m]
)

factorial[6]

Another way to program this function with Recursive calls , or calling the same program inside the program.（ 再帰呼び出し)

factorialR[n_] := (
  If[n == 1  ,  Return[n]  ,  Return[  factorialR[n - 1]*n]   ];
  )

factorialR[6]

Exercice:

Write a program that takes as input k and n  and returns the binomial coefficient  nCk = n! / k!(n-k)!

binom[ n_ , k_ ] :=(

Return[  n!/((k!(n-k)!))];
)

binom[5,3]

Binomial[5,3]

Newton’s method　ニュートン法

Refer to L4-3 page 17 for the definition of Newton’s method to approximate a solution of f(x)=0 
f(x)=0の解を近似するためのニュートン法の定義するには、L4-3のページ17 を参照してください。

Input (入力）: 	differentiable function f   (微分可能な関数　f)
	starting value x0  (初期値　x0)   
	number of iteration  “iter”  (反復の個数  iter )
Output　(出力）: approximation of a solution

newton[f_,x0_,iter_ ]:=(
g[t_]=f'[t];   (*  g is the derivative function of f*)
y=x0;   
Do[ z=y-f[y]/g[y]  (*   Newton iteration  *) ; 
 y=z  (*  swith z and y  *) ;  
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
Return[y];)

Test the program:

f[x_]:=Log[x]-1   (* solutin is x=e*)

Choose x0=2 and  try 1 iterations

newton[f,2,1]

Hmmm,  this is exact value....we need numerical approximation -> put N[...] wherever useful.

newton[f_,x0_,iter_ ]:=(
g[t_]=f'[t];   (*  g is the derivative function of f*)
y=N[x0];   
Do[ z=N[y-f[y]/g[y]]  (*   Newton iteration  *) ; 
 y=N[z]  (*  swith z and y  *) ;  
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
Return[y];)

newton[f,2,1]

Much better! let’s try 2, 3 ,4 iterations

newton[f,2,2]
newton[f,2,3]
newton[f,2,4]

Hmmm...it is better but the precision is not sufficient (see question below).

Application to the approximation of  (Sqrt[2])

Question 1: Find a function g[x_] which has (Sqrt[2])for solution ?
 解(Sqrt[2])を持つ関数ｇを求めよ。

f[x_]:=x^2-2

newton[f,1.4,2]  (* starting value equal to x0=1.4 ,   number of iterations equal to 2 *)
 newton[f,1.4,3(* number of iterations: 3*) ]

Mathematica approximates to 4 digits by default. To modify it for example π up to 32 digits, we can use :  N[ Pi  , 32 ]
Mathematicaは近似するディフォルト値を４桁とする。この値を変更するのに、例えばここで３２桁までのπの近似、　次のコマンドを使う： N[Pi,32]

N [ Pi ,32]

N[Pi,100]

s2=N[Sqrt[2] , 32]

Question 2: modify the program newton[f_ , x0_ , iter_]  as follows.
上のプログラム newton[f_ , x0_ , iter_]　を次のように変更せよ。
Add a 4th input parameter for approximation  “approx_”. And use N [  .... . approx]  where necessary.
明確さを選べるように、４番の入力パラメーター “approx_” を加えよ。それから、必要に応じて、N[ ... , approx]を記入する。

newton[f_,x0_,iter_, approx_ ]:=(
g[t_]=f'[t];   (*  g is the derivative function of f*)
y=N[x0,approx];   
Do[ z=N[y-(f[y]/g[y]) , approx]  (*   Newton iteration  *) ; 
 y=N[z , approx]  (*  switch z and y  *) ;  
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
Return[y];)

newton[f,1.4,3,32]   (* put precision to 32 digits *)

Hmm... still not 32 digit precision ! Why ?? Thisi because the starting value x0=1.4  is input at 5 digits. Let us try 14/10 instead

newton[f,14/10,3,32]

We finally get thesee 32 digits.

Precision: 精度

If two successive approximations are close up to 10E-6 then it is reasonnable to think that the last one has 6 digits precision
もしも相次いでいる二つ近似が10E-6程度近ければ、もっともに6桁が正しいと考えてもいい。

newton[f_,x0_,iter_, approx_ ,prec_]:=(
g[t_]=f'[t];   (*  g is the derivative function of f*)
y=N[x0,approx]; 
stop=False;
Do[ z=N[y-(f[y]/g[y]) , approx] (*   Newton iteration  *) ; Print[z];
If[Abs[y-z]<prec,Return[stop=True]];
 y=N[z , approx] (*  switch z and y  *) ;  
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
If[stop,Return[y],Print["Fail to reach precision ",prec]];)

Question 3: Find a approximation of (Sqrt[2]) at 100 digits precision with the program newton (find the number of iterations necesary iter_)
プログラムnewtonを使って(Sqrt[2])の100桁までの近似を求めよ。

f[x_]=x^2-2

newton[f,1,7,110,10^(-22)]

s2=N[Sqrt[2],100]

Slow convergence at a multiple zero.重複根における遅い収束

f[x_]:=(x^2-2)^2

newton[f,1,20,20,10^(-10)]

Back to Problem

We can print the value of the difference between s2=N [Sqrt[2] , 100]  found after each iterations

newton[f_,x0_,iter_, prec_ ]:=(
g[t_]=f'[t];   (*  g is the derivative function of f*)
y=N[x0,prec];   
Do[ z=N[y-(f[y]/g[y]) , prec]  (*   Newton iteration  *) ; 
 y=N[z , prec]  (*  switch z and y  *) ;  
Print[s2-y];
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
Return[y];)

newton[f,14/10,5,100]

After 5 iterations, we find an approximation at 74 digits. Let’s try 6 iterations.

newton[f,14/10,6,100]

Question 4: Consider the  f[x_] := (x^7) - 3 (x^6) - (x^5)^+ 4 (x^4) - 3 (x^3)+2(x^2)-x+1.
Plot it over the range {x,-1,3}  and the PlotRange->{-1,1}.
Find an approximation of the root close to 3.

Secant’s Method 割線法

secant[f_,x0_,x1_,iter_, approx_ ,prec_]:=(
y0=N[x0,approx];
y1=N[x1,approx];  
stop=False;
Do[ 
y2=N[y1-(f[y1](y1-y0)/(f[y1]-f[y0])) , approx]  (*   Newton iteration  *) ; Print[y2]; 
If[
 Abs[y2-y1]<prec , 
 Return[stop =True],
 y0=N[y1,approx];  
 y1=N[y2,approx]  ;
];
, {i,1,iter} (* repeat from 1 to iter *)
 ];
If[stop,Return[y2],Print["Fail after ",iter," iterations. Last two values computed: ",y2," and ",y1,". Difference: ",y2-y1]];)

func[x_]:=x^2-2

secant[func,1,2,20,110,10^(-100)]

a;

Two Variables Newton’s Method 2変数ニュートン法

newton2[F_,X0_,iter_, approx_,prec_ ]:=(
jacobi[xx_,yy_]=Outer[D,F[xx,yy],{xx,yy}]  (*  jac is the derivative function of f*);
Y=Table[N[i,approx],{i,X0}];  
stop=False;
Do[ 
	Z=Y-Inverse[jacobi[Y[[1]],Y[[2]]]].F[Y[[1]],Y[[2]]]  (*   Newton iteration  *) ;
Print[Z];
If[ Norm[Z-Y]<prec,Return[stop=True],
	Y=Table[N[i , approx],{i,Z}]  (*  switch z and y  *) ;  
]
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
If[stop,Return[Z],Print["Fail after ",iter," iterations" ]]
;)

F[x_,y_]={x+y^2-1,x*y+y};

F[xx,yy][[1]]

F/.{x,y}->{1,1}

D[F[xx,yy][[1]],yy]

{Outer[D,F[x,y],{x,y}], D[F[x,y],{y,2}]}

newton2[F,{10,9},10,30,10^(-5)]

Newton in n variables

rule[A_,B_]:=Map[Rule@@#&,Table[{A[[i]],B[[i]]},{i,1,Length[A]}]];

Algorithm with choice of variables (maybe unnecessary for numerical solution)

newtonMV[F_,X0_,var_,iter_, approx_,prec_ ]:=(
jac=Outer[D,F@@var,var] ; (*  jac is the derivative function of f*);
Y=Table[N[i,approx],{i,X0}];  
stop=False;
Do[ 
	Z=Y-Inverse[jac/.rule[var,Y]].F@@Y  (*   Newton iteration  *) ;
Print[Z];
If[ Norm[Z-Y]<prec,Return[stop=True],
	Y=Table[N[i , approx],{i,Z}]  (*  switch z and y  *) ;  
]
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
If[stop,Return[Z],Print["Fail after ",iter," iterations" ]]
;)

newtonMV[F,{10,9},{xx,yy},10,30,10^(-5)]

Algorithm without variables in argument (maybe better in  the numerical case)

newtonM[F_,X0_,iter_, approx_,prec_ ]:=(
n:=Length[X0];
var=Table[x[i],{i,1,n}]; (*only change compared to the above*)
jac=Outer[D,F@@var,var] ; (*  jac is the derivative function of f*);
Y=Table[N[i,approx],{i,X0}];  
stop=False;
Do[ 
	Z=Y-Inverse[jac/.rule[var,Y]].F@@Y  (*   Newton iteration  *) ;
Print[Z];
If[ Norm[Z-Y]<prec,Return[stop=True],
	Y=Table[N[i , approx],{i,Z}]  (*  switch z and y  *) ;  
]
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
If[stop,Return[Z],Print["Fail after ",iter," iterations" ]]
;)

newtonM[F,{10,9},10,30,10^(-5)]

g[x1_,x2_,x3_]={3 x1-Cos[x2 x3]-1/2,x1^2-81(x2+0.1)^2+Sin[x3]+1.06,Exp[-x1 x2]+20x3+1/3(10 Pi-3)}

newtonM[g,{0.1,0.1,-0.1},5,20,10^(-10)]

Implicit Newton Method

rule[A_,B_]:=Map[Rule@@#&,Table[{A[[i]],B[[i]]},{i,1,Length[A]}]];

ImplicitNewton1[F_,var_,Y0_,iter_]:=(
jac1=D[F[var,yy],yy]  ;(*  jac is the derivative function of f*);Print[jac1];
Y=Y0;(*t=var; *)
Do[ 
	Z=Y-(jac1/.yy->Y)^(-1)*F[var,Y]  ; (*   Newton iteration  *) Print[Z];
Z=Normal[Z+O[var]^(2^i)];
Y=Z;
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
Return[Z];)

Clear[ImplicitNewton1]

ImplicitNewton1[f,t,1,2]

f[x_,y_]:=Sin[x*y]+x^2*Cos[Log[1+x]-y]

ImplicitNewton1[f,t,{1},2]

Normal[Series[f[a,b],{a,0,3},{b,0,3}]]

ImplicitNewton[F_,varSeq_,Y0_,p_Integer,iter_]:=(   (* Y0 is a point such that jac_X0(F,var) is invertible -- F must be a sequence !!*)
n:=Length[Y0];
var=Table[x[i],{i,1,n+p}];
varY=Table[var[[i]],{i,varSeq}]; (* variables to which the Jacobean must be taken *)
varX=Complement[var,varY];Print[varX];
jac1=Outer[D,F@@var,varY]  ;Print[jac1//MatrixForm];(*  jac is the derivative function of f*);
Y=Y0;(*t=var; *)
Do[ 
	Z=Y-Inverse[(jac1/.rule[varY,Y])].F@@var/.rule[varY,Y]  ; (*   Newton iteration  *) Print[Z];
Do[
tmp=(Z[[j]]/.rule[varX,t*varX] )+O[t]^(2^i);Print[tmp];
Z[[j]]=Normal[tmp]/.t->1;
(* Z=Table[Normal[(Z[[j]]/.rule[varY,t*varY] )+O[t]^(2^i)]/.t->1,{j,1,n}]; 
Print[Z];*)
,{j,1,n}];
Y=Z;
   , {i,1,iter} (* repeat from 1 to iter *)
 ];
Return[Z];)

F[x_,y_,z_]={x^2+y^2+z^2-1,x^2-y^2+z}

F[a,b,c]

Solve[F[x,y,z]=={0,0},{x,y}]

ImplicitNewton[F,{1,2},{Sqrt[2]/2,Sqrt[2]/2},1,3]

f1[x_,y_,z_]={x^2+y^2+z^2-1}

ImplicitNewton[f1,{1},{1/3},2,3]