Jeremy Knight Math 639 Algorithm 10.2
<Text-field style="Heading 1" layout="Heading 1">gauss(a,b,n)</Text-field> Solve the linear system ax = b where a is an array(1..n,1..n) b is an array(1..n) x is an array(1..n) nrow is the pivoting array(1..n) n = number of variables gauss:=proc(a, b, n) local i, j, k, p, m, tmp, maxval, tmpval, s, x, nrow; x:=Array(1..n); 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; x; end:
<Text-field style="Heading 1" layout="Heading 1">MxInverse(n,A)</Text-field> MxInverse:=proc(n,A) #This algorithm computes the inverse of an nxn matrix using #the gauss algorithm to solve the system Ax=b local i,j,b,x,k,m,Inv,A0: #option trace: A0:=Array(1..n,1..n): Inv:=Array(1..n,1..n): for i from 1 to n do for k from 1 to n do for m from 1 to n do A0[k,m]:=A[k,m]: end do: end do: b:=Array(1..n): b[i]:=1: x:=gauss(A0,b,n): for j from 1 to n do Inv[j,i]:=x[j]: end do: end do: return (Inv): end proc:
<Text-field style="Heading 1" layout="Heading 1">Matrix Vector Product</Text-field> This algorithm will return the product of a matrix and a vector MVProduct:=proc(n,A,x) #This algorithm will return the product of a matrix and a vector local i,j,p; #option trace: p:=Array(1..n); for i from 1 to n do for j from 1 to n do p[i]:=p[i]+A[i,j]*x[j]: end do: end do: return(p): end proc:
JSFH
<Text-field style="Heading 1" layout="Heading 1">Broyden(n,x0,f,Jac,TOL,N)</Text-field> JSFH Broyden:=proc(n,x0,f,Jac,TOL,N) #n=number of equations and unknowns, x0= initial approximation #f=vector of functions, Jac=Jacobian derivative matrix #TOL=tolerance, N=maximum iterations local i,j,k,p,u,x,s,y,v,w,z,temp1,temp2,A,A0,A1,InfNorm: #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: #####end subroutine section##### #define x as the iterate values of the functions and #define J as the Jacobian matrix approximate, and x:=x0; ########J:=Array(1..n,1..n): A0:=Array(1..n,1..n): #####end definitions##### #Set A0=J(x_0) where J(x)=Jacobian; Set v=F(x_0) v:=Array(1..n): for i from 1 to n do for j from 1 to n do A0[i,j]:=Jac[i,j](x0): end do: v[i]:=f[i](x0): end do: ##############print("***********A0=",A0): #Set A1=A0^(-1) using Gaussian elimination to solve Ax=b A1:=MxInverse(n,A0): #############print("********A1,v",A1,v): ###################s:=MVProduct(n,A1,v): s:=(-1)*MVProduct(n,A1,v): x:=x+s: k:=1: print(k,x): k:=2: while k<=N do w:=Array(1..n): for i from 1 to n do w[i]:=v[i]: #Save v end do: for i from 1 to n do v[i]:=f[i](x): #v=F(x^k) end do: y:=v-w: #y=y[k] ############################z:=MVProduct(n,A1,y): z:=(-1)*MVProduct(n,A1,y): temp1:=0: ###################################### for i from 1 to n do temp1:=temp1+s[i]*z[i]: end do: p:=(-1)*temp1: #Set p=s[k]^t*A[k=1]^(-1)*y[k] u:=Array(1..n): for i from 1 to n do: for j from 1 to n do u[i]:=u[i]+s[j]*A1[j,i] end do: end do: temp1:=s+z: temp2:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do temp2[i,j]:=temp1[i]*u[j]: end do: end do: A1:=A1+(1/p)*temp2; #A=A[k]^(-1) #################3print("********** A1=",A1): s:=(-1)*MVProduct(n,A1,v): #s=-A[k]^(-1)F(x[k]) x:=x+s: #x=x[k+1] if InfNorm(n,s)<TOL then print("x=",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:
JSFH f1:=(x::Array)->3*x[1]-cos(x[2]*x[3])-(1.0/2); f2:=(x::Array)->x[1]^2-81*(x[2]+0.1)^2+sin(x[3])+1.06; f3:=(x::Array)->evalf(exp(-x[1]*x[2])+20*x[3]+(10.0*(Pi)-3.0)/3.0); F:=Array([f1,f2,f3]); Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCgmOSQ2IyIiIiIiJC1JJGNvc0c2JEYoSShfc3lzbGliR0YmNiMqJiZGLjYjIiIjRjAmRi42I0YxRjAhIiIqJiQiIzVGPUYwRjpGPUY9RiZGJkYm Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCoqJCY5JDYjIiIiIiIjRjEqJCwmJkYvNiNGMkYxJEYxISIiRjFGMiEjIiktSSRzaW5HNiRGKEkoX3N5c2xpYkdGJjYjJkYvNiMiIiRGMSQiJDEiISIjRjFGJkYmRiY= Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLUkmZXZhbGZHRig2IywoLUkkZXhwRzYkRihJKF9zeXNsaWJHRiY2IywkKiYmOSQ2IyIiIkY6JkY4NiMiIiNGOiEiIkY6JkY4NiMiIiQiIz8qJiwmKiYkIiQrIkY+RjpJI1BpR0YoRjpGOiQiI0lGPkY+RjpGSUY+RjpGJkYmRiY= LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKmtrNWoj J:=Array(1..3,1..3): J[1,1]:=(x::Array)->3: J[1,2]:=(x::Array)->x[3]*sin(x[2]*x[3]): J[1,3]:=(x::Array)->x[2]*sin(x[2]*x[3]): J[2,1]:=(x::Array)->2*x[1]: J[2,2]:=(x::Array)->-162*(x[2]+0.1): J[2,3]:=(x::Array)->cos(x[3]): J[3,1]:=(x::Array)->-x[2]*exp(-x[1]*x[2]): J[3,2]:=(x::Array)->-x[1]*exp(-x[1]*x[2]): J[3,3]:=(x::Array)->20: J; LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKltQcmsj J[3,2](Array([1.0,2.0,3])); JCErS0dOYDghIzU= x0:=Array([0.1,0.1,-0.1]); Broyden(3,x0,F,J,10^(-10),20); LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKi82QWsj NiQiIiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqK3EwaSM= NiQiIiMtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqWzhYXiM= NiQiIiQtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqc2FgayM= NiQiIiUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqZ0woZkU= NiQiIiYtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqNzNpaSM= NiQiIictSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqRzF4XiM= NiQiIigtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqSzVfXiM= NiQiIiktSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqP2psaiM= NiRRI3g9NiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEdGJCIqU1k2ayM= NiRRLGl0ZXJhdGlvbnM9NiIiIio= JSFH JSFH
<Text-field style="Heading 1" layout="Heading 1">Broyden2(n,x0,f,TOL,N)</Text-field> LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= Broyden2:=proc(n,x0,f,TOL,N) #n=number of equations and unknowns, x0= initial approximation #f=vector of functions, Jac=Jacobian derivative matrix #TOL=tolerance, N=maximum iterations local i,j,k,p,u,x,s,y,v,w,z,temp1,temp2,A,A0,A1,InfNorm,e,h: #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: #####end subroutine section##### #define x as the iterate values of the functions and #define J as the Jacobian matrix approximate, and x:=x0; ########J:=Array(1..n,1..n): A0:=Array(1..n,1..n): #####end definitions##### # Set v=F(x_0) v:=Array(1..n): for i from 1 to n do v[i]:=f[i](x0): end do: #calculate the Jacobian J(x). Set A0=J(x_0) where J(x)=Jacobian; e:=Array(1..n): h:=.001: for i from 1 to n do for j from 1 to n do e[j]:=h: A0[i,j]:=(f[i](x+e)-f[i](x-e))/(2*h): e[j]:=0: end do: end do: ##############print("***********A0=",A0): #Set A1=A0^(-1) using Gaussian elimination to solve Ax=b A1:=MxInverse(n,A0): #############print("********A1,v",A1,v): ###################s:=MVProduct(n,A1,v): s:=(-1)*MVProduct(n,A1,v): x:=x+s: k:=1: print(k,x): k:=2: while k<=N do w:=Array(1..n): for i from 1 to n do w[i]:=v[i]: #Save v end do: for i from 1 to n do v[i]:=f[i](x): #v=F(x^k) end do: y:=v-w: #y=y[k] ############################z:=MVProduct(n,A1,y): z:=(-1)*MVProduct(n,A1,y): temp1:=0: ###################################### for i from 1 to n do temp1:=temp1+s[i]*z[i]: end do: p:=(-1)*temp1: #Set p=s[k]^t*A[k=1]^(-1)*y[k] u:=Array(1..n): for i from 1 to n do: for j from 1 to n do u[i]:=u[i]+s[j]*A1[j,i] end do: end do: temp1:=s+z: temp2:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do temp2[i,j]:=temp1[i]*u[j]: end do: end do: A1:=A1+(1/p)*temp2; #A=A[k]^(-1) #################3print("********** A1=",A1): s:=(-1)*MVProduct(n,A1,v): #s=-A[k]^(-1)F(x[k]) x:=x+s: #x=x[k+1] if InfNorm(n,s)<TOL then print("x=",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:
JSFH
<Text-field style="Heading 1" layout="Heading 1">Exercise set 10.3</Text-field>
<Text-field style="Heading 2" layout="Heading 2">1a</Text-field> JSFH f1:=(x::Array)->4*x[1]^2-20*x[1]+1/4*x[2]^2+8; f2:=(x::Array)->1/2*x[1]*x[2]^2+2*x[1]-5*x[2]+8; F:=Array([f1,f2]); x0:=Array([0,0]); Broyden2(2,x0,F,.001,10); Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCoqJCY5JDYjIiIiIiIjIiIlRi4hIz8qJCZGLzYjRjJGMiNGMUYzIiIpRjFGJkYmRiY= Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCoqJiY5JDYjIiIiRjEmRi82IyIiI0Y0I0YxRjRGLkY0RjIhIiYiIilGMUYmRiZGJg== LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKiFbMFhF LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKldbJUdF NiQiIiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqdzlLayM= NiQiIiMtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqWyMqel0j NiQiIiQtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqSyQpb2sj NiRRI3g9NiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEdGJCIqJVs1UkU= NiRRLGl0ZXJhdGlvbnM9NiIiIiU=
JSFH TTdSMApJNlJUQUJMRV9TQVZFLzI2MzEwNjQ2NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCUjZjFHJSNmMkclI2YzR0YmTTdSMApJNlJUQUJMRV9TQVZFLzI2NDcxMzc0OFgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIyoiJCIkZio2IyclInhHJSZBcnJheUdGJjYkJSlvcGVyYXRvckclJmFycm93R0YmIiIkRiZGJkYmZipGKEYmRixGJiwkJjkkNiMiIiIiIiNGJkYmRiZmKkYoRiZGLEYmLCQqJiZGMzYjRjZGNS0lJGV4cEc2IywkKiZGMkY1RjpGNSEiIkY1RkFGJkYmRiZmKkYoRiZGLEYmKiYmRjM2I0YvRjUtJSRzaW5HNiMqJkY6RjVGREY1RjVGJkYmRiZmKkYoRiZGLEYmLCZGOiEkaSIkRjVGQUZMRiZGJkYmZipGKEYmRixGJiwkKiZGMkY1RjxGNUZBRiZGJkYmZipGKEYmRixGJiomRjpGNUZGRjVGJkYmRiZmKkYoRiZGLEYmLSUkY29zRzYjRkRGJkYmRiZmKkYoRiZGLEYmIiM/RiZGJkYmRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2NDIyMTEwNFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiIiIhIiJGJyRGKUYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2MjA1NzAwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK0RucCkqXCEjNSQiK19bb1k+ISM2JCErPlo/Ol9GKUYmTTdSMApJNlJUQUJMRV9TQVZFLzI1MTQ1MTM0OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK2BQJykqKlwhIzUkIipJUnl0KSEjNiQhK1hkdUpfRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2NDUzNTQ3MlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK3JmMStdISM1JCIqek5GbikhIzckISs5TXNOX0YpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2NTk3MzM2MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKylHLisrJiEjNSQiKiRSJEcmUiEjOCQhK15vKGZCJkYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2MjYyMDgxMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKzwrKytdISM1JCIpMnpNPiEjOSQhKy54KWZCJkYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTc3MDYyOFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKytdISM1JCIoXUxnIiEjOyQhK154KWZCJkYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTUyMTAzMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKytdISM1JCErSVpNdTghIz4kIStjeClmQiZGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzI2MzY1NjMyMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKytdISM1JCEqIlJjJGYjISM+JCErY3gpZkImRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2NDExNDY0MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKytdISM1JCEqIlJjJGYjISM+JCErY3gpZkImRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2NDUwNTQ4MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiMiIyUjZjFHJSNmMkdGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2Mjg0NDg0NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiMiIyIiIUYnRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2NDMyMTQ3NlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiMiIyQiKysrKytTISM1JCIrKysrZzwhIipGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI1MDc5OTI0OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiMiIyQiK0EsI3p4JSEjNSQiK0U3VEY+ISIqRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2NDY4ODMzMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiMiIyQiKzlPODRdISM1JCIrRmoiKioqPiEiKkYmTTdSMApJNlJUQUJMRV9TQVZFLzI2MzkxMDQ4NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiMiIyQiK0A0Ny1dISM1JCIrPFdfKz8hIipGJg==