As an example, the common intersection point of 3 spheres in 3-dim space is calculated.
The 3 spheres each have an equal radius of 7, and their midpoints (x,y,z) are (3,3,3), (1,-2,-1) and (-1,2,-2) respectively. These spheres have a common intersection point. If you enlarge the radius of one of the spheres such that it encompasses the other two, the result "no solution" will be given.
The problem-related data is in yellow, the rest of the program remains unchanged for any problem.
Code: Select all
option base 1
dim jac(9,9),dx(9),eps(9),var(9)
max_iter=10 ! max_err=.01
'y' ' # of variables and starting values for them
n=3 ! read var(1),var(2),var(3) ! data 0,0,0
''
for iter=1 to max_iter
for i=1 to n
eps(i)=func(i,var) ' deviations from zero to be corrected
for j=1 to n ! jac(i,j)=deriv(i,j,var) ! next j ' Jacobian matrix
next i
lin_eqn(n,jac,dx,eps) ! err=0 ' corrections to the variables
for i=1 to n ! var(i)-=dx(i) ! err+=dx(i)^2 ! next i ' apply them
if err<max_err then break
next iter
if err>max_err then ! print "no solution"
else ! for i=1 to n ! print "var"&i&" : ";var(i) ! next i
end if
end
def func(i,x())
'y'
on i goto sphere1,sphere2,sphere3
sphere1: return (x(1)-3)^2+(x(2)-3)^2+(x(3)-3)^2-7^2
sphere2: return (x(1)-1)^2+(x(2)+2)^2+(x(3)+1)^2-7^2
sphere3: return (x(1)+1)^2+(x(2)-2)^2+(x(3)+2)^2-7^2
''
end def
def deriv(i,j,x()) ' partial derivate of function i w.r.t. variable j
fx=func(i,x)
delta=.001*x(j) ! delta=max(abs(delta),.001)
x(j)+=delta ! fxplus=func(i,x) ! x(j)-=delta
return (fxplus-fx)/delta
end def
def lin_eqn(nvar,a(,),x(),b()) ' solve nvar lineair equations
for i=1 to nvar-1
for j=i+1 to nvar
fac=a(j,i)/a(i,i) ! b(j)=b(j)-fac*b(i)
for k=i+1 to nvar ! a(j,k)=a(j,k)-fac*a(i,k) ! next k
next j
next i
x(nvar)=b(nvar)/a(nvar,nvar)
for i=nvar-1 to 1 step -1 ! x(i)=b(i)
for j=i+1 to nvar ! x(i)=x(i)-a(i,j)*x(j) ! next j
x(i)=x(i)/a(i,i)
next i
return
end def