Jeremy Knight Math 639 3/11/11 Algorithm 10.1 Newton's Method for systems To approximate the solution of the nonlinear system LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2J1EiRkYnLyUlYm9sZEdRJXRydWVGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJWJvbGRGJy8lK2ZvbnR3ZWlnaHRHRjctSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RicvRjNGMS9GNlEnaXRhbGljRicvRjZRJ25vcm1hbEYnRkUtSSNtb0dGJDYtUSI9RidGRS8lJmZlbmNlR0Y0LyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRmVuLUkjbW5HRiQ2JlEiMEYnRi9GNUY4RkU= given an initial approximation x. Inputs: newton(f, n, x0, tol) where f is a vector of functions returning a vector of values, J is a matrix of functions which returns a matrix of values, n is the number of equations, x0 is the starting approximation, and tol is the desired tolerance. We will compute an approximation for the Jacobian matrix LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2J1EiSkYnLyUlYm9sZEdRJXRydWVGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJWJvbGRGJy8lK2ZvbnR3ZWlnaHRHRjctSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RicvRjNGMS9GNlEnaXRhbGljRicvRjZRJ25vcm1hbEYnRkVGRQ== by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiSkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRI2lqRidGMkY1RjJGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEieEYnRjJGNS9GNlEnbm9ybWFsRidGSC1JI21vR0YkNi1RIn5GJ0ZILyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZQLyUpc3RyZXRjaHlHRlAvJSpzeW1tZXRyaWNHRlAvJShsYXJnZW9wR0ZQLyUubW92YWJsZWxpbWl0c0dGUC8lJ2FjY2VudEdGUC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmluLUZLNi1RKSZhcHByb3g7RidGSEZORlFGU0ZVRldGWUZlbi9GaG5RLDAuMjc3Nzc3OGVtRicvRltvRmBvLUkmbWZyYWNHRiQ2KC1GQTYkLUYjNiotRiw2JS1GLzYlUSJGRidGMkY1LUYjNiUtRi82JVEiaUYnRjJGNUYyRjVGPS1GQTYkLUYjNihGRS1GSzYtUSIrRidGSEZORlFGU0ZVRldGWUZlbi9GaG5RLDAuMjIyMjIyMmVtRicvRltvRltxLUYvNiVRImhGJ0YyRjUtRks2LVEnJnNkb3Q7RidGSEZORlFGU0ZVRldGWUZlbkZnbkZqbi1GLDYlLUYvNiVRImVGJ0YyRjUtRiM2JS1GLzYlUSJqRidGMkY1RjJGNUY9RkhGSEZKLUZLNi1RKiZ1bWludXMwO0YnRkhGTkZRRlNGVUZXRllGZW5GanBGXHFGSkZpby1GQTYkLUYjNihGRS1GSzYtUSgmbWludXM7RidGSEZORlFGU0ZVRldGWUZlbkZqcEZccUZdcUZgcUZjcUZIRkhGSEZILUYjNictSSNtbkdGJDYkUSIyRidGSEZgcUZdcUYyRjUvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmJzLyUpYmV2ZWxsZWRHRlBGSA==.
<Text-field style="Heading 1" layout="Heading 1">newton(f,n,x0,tol,N)</Text-field> LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= newton:=proc(n,f,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: ##########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: #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:
g1:=(x::Array) -> -cos(x[2]*x[3])-(1/2); f1:=(x::Array)->g1(x)+3.0*x[1]; f2:=(x::Array)->x[1]^2-81.0*(x[2]+0.1)^2+sin(x[3])+1.06; f3:=(x::Array)->evalf(exp(-1.0*x[1]*x[2])+20.0*x[3]+(10.0*Pi-3.0)/3.0); f:=Array([f1,f2,f3]); Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCYtSSRjb3NHNiRGKEkoX3N5c2xpYkdGJjYjKiYmOSQ2IyIiIyIiIiZGNDYjIiIkRjchIiIjRjtGNkY3RiZGJkYm Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCYtSSNnMUdGJjYjOSQiIiIqJiQiI0khIiJGMSZGMDYjRjFGMUYxRiZGJkYm Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCoqJCY5JDYjIiIiIiIjRjEqJiQiJDUpISIiRjEqJCwmJkYvNiNGMkYxJEYxRjZGMUYyRjFGNi1JJHNpbkc2JEYoSShfc3lzbGliR0YmNiMmRi82IyIiJEYxJCIkMSIhIiNGMUYmRiZGJg== Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLUkmZXZhbGZHRig2IywoLUkkZXhwRzYkRihJKF9zeXNsaWJHRiY2IywkKigkIiM1ISIiIiIiJjkkNiNGOkY6JkY8NiMiIiNGOkY5RjoqJiQiJCsjRjlGOiZGPDYjIiIkRjpGOiomLCYqJiQiJCsiRjlGOkkjUGlHRihGOkY6JCIjSUY5RjlGOkZNRjlGOkYmRiZGJg== LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKmckXDdF x0:=Array([0.1,0.1,-0.1]); newton(3,f,x0,.0000001,20); LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKks/RGgj NiQiIiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqS2pNaCM= NiQiIiMtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqJSk+V2gj NiQiIiQtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqQzRhaCM= NiQiIiUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqWz1raCM= NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKiM+UTxF NiRRLGl0ZXJhdGlvbnM9NiIiIiY= x0:=Array([0.1,0.1,-0.1]); newton(3,f,x0,.0000001,20); LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKmtka1wj NiQiIiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqJVtEKVwj NiQiIiMtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqc2MlKlwj NiQiIiQtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqU00xXSM= NiQiIiUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqZ2g9XSM= NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKnchKUhdIw== NiRRLGl0ZXJhdGlvbnM9NiIiIiY=
<Text-field style="Heading 1" layout="Heading 1">newton2(n,f,Jac,x0,tol,N)</Text-field> JSFH JSFH JSFH newton2:=proc(n,f,Jac,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: #####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): 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 J[i,j]:=Jac[i,j](x): end do: end do: #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: