Jeremy Knight Math 639 Newtons Minimum Method
<Text-field style="Heading 1" layout="Heading 1">NewtonMin(f,n,x0,tol,N)</Text-field> JSFH NewtonMin:=proc(n,f,delf,x0,TOL,N) #f = vector of functions, n = number of equations and variables, #x0=initial approximation, tol=tolerance #f must be of the form >>f=(x::array)->function_definition #M= maximum number of iterations local e,h,F,J,i,j,k,x,y,InfNorm,gauss,temp: #option trace; #####Subroutine section##### ##InfNorm for finding norm for tolerance InfNorm:=proc(n,x) local inf_norm,j,p: inf_norm:=x[1]: for j from 2 to n do if abs(x[j])>inf_norm then inf_norm:=abs(x[j]): p:=j: end if: end do: return(inf_norm); end proc: ##gauss - for solving J(x)*y=-F(x) gauss:=proc(A, bb, n) local a,b,i, j, k, p, m, tmp, maxval, tmpval, s, x, nrow; x:=Array(1..n); #Copy A and bb to a and b to prevent original arrays from being modified a:=Array(1..n,1..n): b:=Array(1..n): for i from 1 to n do for j from 1 to n do a[i,j]:=A[i,j]: end do: b[i]:=bb[i]: end do: nrow:=Array(1..n); for i from 1 to n do nrow[i]:= i; end do; s:=Array(1..n); for i from 1 to n do maxval:=abs(a[i,1]); for j from 2 to n do tmpval:=abs(a[i,j]); if ( tmpval > maxval ) then maxval:=tmpval; end if; end do; s[i]:=maxval; end do; for i from 1 to n-1 do # find pivot row maxval:= abs(a[nrow[i],i])/s[nrow[i]]; p:=i; for j from i+1 to n do tmpval:= abs(a[nrow[j],i])/s[nrow[j]]; if (tmpval > maxval) then p:=j; maxval:=tmpval; end if; end do; if ( maxval = 0 ) then return(0); end if; # switch rows if ( p <> i ) then tmp:=nrow[p]; nrow[p]:=nrow[i]; nrow[i]:=tmp; end if; # zero out column i for k from i+1 to n do m:=a[nrow[k],i]/a[nrow[i],i]; a[nrow[k],i]:=0; for j from i+1 to n do a[nrow[k],j]:=a[nrow[k],j]-m*a[nrow[i],j]; end do; b[nrow[k]]:=b[nrow[k]]-m*b[nrow[i]]; end do; if ( a[nrow[n],n] = 0 ) then return(0); end if; end do; # backward substitution x[n]:= b[nrow[n]]/a[nrow[n],n]; for j from n-1 to 1 by -1 do x[j]:=b[nrow[j]]; for k from n to j+1 by -1 do x[j]:=x[j]- a[nrow[j],k]*x[k]; end do; x[j]:=x[j]/a[nrow[j],j]; end do; return(x); end: gradient:=proc(x::Array) local i,j,F,J,temp,DEL: #del{g(x)}=2*J(x)^t*F(x); #First calculate the Jacobian J(x) #J:=Jacobian(x): DEL:=Array(1..n): #for i from 1 to n do # F[i]:=f[i](x): # end do: for i from 1 to n do DEL[i]:=delf[i](x): end do: return(DEL): end proc: ##########end subroutine section########## #define x as the iterate values of the functions and #define J as the Jacobian matrix approximate, and #e = vector representing the jth column of the idenity matrix #h = step value for the Jacobian approximation x:=Array(1..n): for i from 1 to n do x[i]:=x0[i]: end do: #J:=Array(1..n,1..n): F:=Array(1..n): e:=Array(1..n): h:=.001: k:=1: #####end definitions##### while k<=N do #calculate F(x)=x for i from 1 to n do F[i]:=f[i](x); end do; # #calculate the Jacobian J(x) # for i from 1 to n do # for j from 1 to n do # e[j]:=h: # J[i,j]:=(f[i](x+e)-f[i](x-e))/(2*h): # e[j]:=0: # end do: # end do: #Calculate the gradiant of f G:=gradient(x): #Solve J(x)y=-F(x) y:=gauss(J,((-1)*F),n): x:=x+y: if InfNorm(n,y)<=TOL then print("solution=", x): print("iterations=",k): break: end if: print(k,x); k:=k+1 end do: if k>N then print("maximum number of iterations exceeded"): end if: end proc: