# # Do all the initialisations # theta:=Array([161.2, 45.10, 309.0]): sigma:= Array([0.8, 0.6, 1.3, 2.0]): X:= Array([746, 629, 1571, 155]): Y:= Array([1393, 375, 259, 987]): d4:=864.3: pi:=3.1415926540: for i to 3 do theta[i]:= 2*pi*theta[i]/360; sigma[i]:= 2*pi*sigma[i]/360; end do: # # Generate equations in appropriate form # for solution # S:=add(((arctan(x-X[i],y-Y[i])-theta[i])/sigma[i])^2, i=1..3) +((((x-X)^2+(y-Y)^2)^(1/2)-d4)/sigma)^2: eqns:={diff(S,x)=0, diff(S,y)=0}: # # Attempt solution # sol:=fsolve(eqns, {x, y}); whattype(sol); # # Hmmm no solutions. Note the denominators: whenever y=Y[i] # looks like one ends up with "0*Infinity" - pretty difficult # to work with. # # Plot the surfaces just to see if there are regions where # solutions look possible # plot3d( [ lhs(eqns), lhs(eqns) ], x = -1000..1000, y = -1000..1000 ); # # Zoom in to a likely region - the big spikey bit # plot3d( [ lhs(eqns), lhs(eqns) ], x = 100..400, y = 200..500 ); # # Looks promising so test for a solution in this region # (Probably many other solutions - unfortunately I suspect # most of them are close to the areas where the functions # become indeterminate # sol:=fsolve ( { diff(S,x) = 0, diff(S,y) = 0}, { x = 100..300, y = 300..500} ); # # *Looks* as if all solutions are for x>0, y>0, and probably(?) # don't extend much beyond max(X), max(Y) so split this region # into sub-regions and search within each subregion # t1:=time(): eqns:= { diff(S,x) = 0, diff(S,y) = 0 }: for xr from 0 by 200 to 2000 do for yr from 0 by 200 to 2000 do sol:= fsolve ( eqns, { x = xr..xr+200, y = yr..yr+200} ); if whattype(sol)=set then printf(" sol=%a\134n", sol); fi; od; od; t2:=time()-t1;