Jeremy Knight Math 639 3/11/11 The Implicit Trapezoidal Method for solving a system of differential equations using Newton's Method to solve a system of nonlinear Equations.
<Text-field style="Heading 1" layout="Heading 1">newton(n,f,x0,tol,N,display)</Text-field> LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= newton:=proc(n,f,x0,TOL,N,display) #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 if display=1 then print("solution=", x): print("iterations=",k): end if: return(x): break: end if: if display=1 then print(k,x); end if: k:=k+1 end do: if k>N then print("maximum number of iterations exceeded"): end if: end proc:
<Text-field style="Heading 1" layout="Heading 1">ImplicitTrapezoidal(n,N,f,a,b,alpha,TOL,M,display)</Text-field> ImplicitTrapizoidal:=proc(n,N,f,a,b,alpha,TOL,M,display) #Uses the Implicit Trapezoidal Method solve a system of differential #equations using Newton's Method. #Must be accompanied by newton(f,n,x0,tol,N) #n=number of equations and variables,N=number of subintervals #F=Array of differential equations in terms of t,x[1],x[2],...,x[n] #time interval=[a,b] #alpha= initial values #TOL=tolerance, M=maximum iterations for Newton local i,j,w,w0,w1,g1,g2,g3,G,t,h,jNewton: #define t[i] as time array and divide into N subintervals t:=Array(1..N): h:=(b-a)/N: t:=a: #define w0 as the previous vector [x_j,y_j,z_j] w0:=Array(1..n): for i from 1 to n do w0[i]:=alpha[i]: end do: print(0,t,w0): #define F as a process to calculate the values of f(x) ##F:=proc(x::Array,t) ## local i,temp: ## temp:=Array(1..n): ## for i from 1 to n do ## temp[i]:=f[i](x): ## end do: ## return(temp): ##end proc: for j from 1 to N do t:=a+j*h: #Define system of trapezoidal equations G(t,x,y,z)=G(t,w[1],w[2],w[3]) g1:=(w::Array)->w0[1]+(h/2.0)*(f[1](w)+f[1](w0))-w[1]: g2:=(w::Array)->w0[2]+(h/2.0)*(f[2](w)+f[2](w0))-w[2]: g3:=(w::Array)->w0[3]+(h/2.0)*(f[3](w)+f[3](w0))-w[3]: G:=Array([g1,g2,g3]): #define w1 as the next approximation w1:=newton(n,G,w0,TOL,M,display): print(j,t,w1): #Let w1 be an array of current (j+1) approximations for i from 1 to n do w0[i]:=w1[i]: end do: end do: end proc: f1:=(x::Array)->-0.013*x[1]-1000.0*x[1]*x[3]; f2:=(x::Array)->-2500.0*x[2]*x[3]; f3:=(x::Array)->-0.013*x[1]-1000.0*x[1]*x[3]-2500.0*x[2]*x[3]; F:=Array([f1,f2,f3]); Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCYqJiQiIzghIiQiIiImOSQ2I0YxRjEhIiIqKCQiJisrIkY1RjFGMkYxJkYzNiMiIiRGMUY1RiZGJkYm Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCQqKCQiJitdIyEiIiIiIiY5JDYjIiIjRjEmRjM2IyIiJEYxRjBGJkYmRiY= Zio2IydJInhHNiJJJkFycmF5RyUqcHJvdGVjdGVkR0YmNiRJKW9wZXJhdG9yR0YmSSZhcnJvd0dGJkYmLCgqJiQiIzghIiQiIiImOSQ2I0YxRjEhIiIqKCQiJisrIkY1RjFGMkYxJkYzNiMiIiRGMUY1KigkIiYrXSNGNUYxJkYzNiMiIiNGMUY5RjFGNUYmRiZGJg== LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKldjamoj ImplicitTrapizoidal(3,5,F,0,50,Array([1,1,0]),.1,15,1); NiUiIiFGIy1JJkFycmF5RyUqcHJvdGVjdGVkRzYjL0kkJWlkRzYiIiojPi9ARA== NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKlsqSEhF NiRRLGl0ZXJhdGlvbnM9NiIiIiI= NiUiIiIiIzUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqWypISEU= NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKj9wU18j NiRRLGl0ZXJhdGlvbnM9NiIiIiI= NiUiIiMiIz8tSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqP3BTXyM= NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKkcmcEdF NiRRLGl0ZXJhdGlvbnM9NiIiIiI= NiUiIiQiI0ktSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqRyZwR0U= NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKiFlXkxF NiRRLGl0ZXJhdGlvbnM9NiIiIiI= NiUiIiUiI1MtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqIWVeTEU= NiRRKnNvbHV0aW9uPTYiLUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHRiQiKm9AJT5E NiRRLGl0ZXJhdGlvbnM9NiIiIiI= NiUiIiYiI10tSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqb0AlPkQ= JCErSCdvbUInISM6
LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= JSFH JSFH LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= TTdSMApJNlJUQUJMRV9TQVZFLzI2MzYzNTY0NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCUjZjFHJSNmMkclI2YzR0YmTTdSMApJNlJUQUJMRV9TQVZFLzI1MjEwNDE5MlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK0JYQiYpZiEjNSQiK0suWiw5ISIqJCErSCdvbUInISM6RiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2MjkyOTk0OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKzAnM0U2KiEjNSQiK1Q/dCkzIiEiKiQhK3lLYyk0KCEjOkYmTTdSMApJNlJUQUJMRV9TQVZFLzI2MjkyOTk0OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKzAnM0U2KiEjNSQiK1Q/dCkzIiEiKiQhK3lLYyk0KCEjOkYmTTdSMApJNlJUQUJMRV9TQVZFLzI1MjQwNjkyMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK2IlKilcQSkhIzUkIis0QV14NiEiKiQiKy92XmE2ISM6RiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MjQwNjkyMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK2IlKilcQSkhIzUkIis0QV14NiEiKiQiKy92XmE2ISM6RiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2Mjg2OTUyOFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyQ+KlJOdSEjNSQiKyJbYGtEIiEiKiQhK2kyImZmJyEjOkYmTTdSMApJNlJUQUJMRV9TQVZFLzI2Mjg2OTUyOFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyQ+KlJOdSEjNSQiKyJbYGtEIiEiKiQhK2kyImZmJyEjOkYmTTdSMApJNlJUQUJMRV9TQVZFLzI2MzM1MTU4MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK3psa2htISM1JCIrQHUkUUwiISIqJCIrW08sdj8hIzpGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MzM1MTU4MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK3psa2htISM1JCIrQHUkUUwiISIqJCIrW08sdj8hIzpGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI1MTk0MjE2OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK0JYQiYpZiEjNSQiK0suWiw5ISIqJCErSCdvbUInISM6RiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTk0MjE2OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK0JYQiYpZiEjNSQiK0suWiw5ISIqJCErSCdvbUInISM6RiY=