MoMath sandbox for program audiences

Project: Math 367-17S
Views: 7
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu1804
# Modeling dropping an object from a height

g = 2 #acceleration of "gravity"
y = 0 #starting y value.  As we move downward, y will INCREASE.  The variable y is POSITION
v = 0 #the velocity starts at zero, since we are just "letting go" of the object
t = 0 #starting time
dt = 0.1 # this is our "delta t" which should be SMALL as we approximate the object's motion and take "snapshots" every dt

# we are going to store the points (t, y) in a list to plot later
data = [(t,y)]  #so the data holds just [(0,0)] to start.

#Now we make a table of t, y, v.
nsteps = 200 #the number of snapshots.  So we will go for nsteps*dt time units.
print("t      y      v")  # the top row

for k in range(nsteps):    #this means that k goes from 0 until nsteps-1 and we loop
#print n(t, digits =3), n(y, digits =3), n(v, digits =3)   #self explanatory (digits tells it how many decimal places to display)
# now we iterate for the next snapshot
t = t + dt            # time advances
y = y + v*dt          # the "delta y" is equal to v*dt, since distance = rate times time!
v = v + g*dt          # same idea
data.append((t,y))    # add the new point to the list of points
list_plot(data)           # when the table is done, plot the t vs. y graph

t y v
# Modeling population growth
#
#  we will assume that P'(t) = aP(t), where a is a smallish constant

a = 0.2 #constant
P = 100 #starting P value.
v = P*a #the starting velocity in people per time units (say, years)
t = 0 #starting time
dt = 0.1 # this is our "delta t" which should be SMALL as we approximate the object's motion and take "snapshots" every dt

# we are going to store the points (t, y) in a list to plot later
data = [(t,P)]  #so the data holds just [(0,0)] to start.

#Now we make a table of t, y, v.
nsteps = 200 #the number of snapshots.  So we will go for nsteps*dt time units.

print("t      P      v")  # the top row

for k in range(nsteps):    #this means that k goes from 0 until nsteps-1 and we loop
#print n(t, digits =3), n(P, digits =3), n(v, digits =3)   #self explanatory (digits tells it how many decimal places to display)
# now we iterate for the next snapshot
t = t + dt            # time advances
P = P + v*dt          # the "delta P" is equal to v*dt, since distance = rate times time!
v = P*a               # same idea
data.append((t,P))    # add the new point to the list of points
list_plot(data)           # when the table is done, plot the t vs. y graph

t P v


#catenary simulation
# we will try to create an approximate solution to the differential equation
# tan(t) = s, where t is theta and s is "ell" (this is a standard letter, and "ell" is hard to read in this font.)

t=0   #start at the bottom of the chain
x=0
y=0
ds =0.01 # this is what we increment length by; a smaller value yields a better approx
s = 0 #the starting arc length from (0,0).  It will increase as we move to the right.
nsteps = 300  #self explanatory

data=[(x,y)]  # So we start with (0,0) and then we will append and finally plot
for k in range(nsteps):   #the expression "range(nsteps)" just means, "the numbers from 0 to nsteps-1, inclusive"
#we will move by a length of ds
# in the direction with slope = tan(t)
# and tan(t)=s, the total distance so far
x += ds*cos(t)     #another programming expression: "x += u" means "increase x by the value u;" i.e. x = x+u
y += ds*sin(t)

#include the printout for debugging
#print k, x, y,t
#append the point
data.append((x,y))

#now increment
s += ds    #computing the total length so far
t = atan(s)  # the angle t is such that tan(t) = s; i.e., the SLOPE equals s at this point

#loop is over, now plot
list_plot(data,aspect_ratio=1)


# compare this with a "real" catenary
#you need to state that x is a variable used in the plot command

%var x
approx = list_plot(data)
real = plot((exp(x)+exp(-x))/2-1, (x,0,1.8), color='red')    #exp(x) means e^x
show(approx+real,aspect_ratio=1)



#square-wheeled trike simulation
# we will try to create an approximate solution to the differential equation
# 1- tan(t) = s, where t is theta and s is "ell" (this is a standard letter, and "ell" is hard to read in this font.)

t=pi/4   #start at the left of the surface, with an angle of 45 degrees (square is on its corner)
x=-.8814 #this is the x-coordinate where the actual curve starts.  We could start with other values if we wanted.
y=0
ds =0.1 # this is what we increment length by; a smaller value yields a better approx
s = 0 #the starting arc length from (0,0).  It will increase as we move to the right.
nsteps = 30  #self explanatory

data=[(x,y)]  # So we start with (0,0) and then we will append and finally plot
for k in range(nsteps):   #the expression "range(nsteps)" just means, "the numbers from 0 to nsteps-1, inclusive"
#we will move by a length of ds
# in the direction with slope = tan(t)
# and tan(t)=s, the total distance so far
x += ds*cos(t)     #another programming expression: "x += u" means "increase x by the value u;" i.e. x = x+u
y += ds*sin(t)

#include the printout for debugging
#print k, x, y,t
#append the point
data.append((x,y))

#now increment
s += ds    #computing the total length so far
t = atan(1-s)  # the angle t is such that 1-tan(t) = s; i.e., the SLOPE equals s at this point

#loop is over, now plot
list_plot(data,aspect_ratio=1)

def spiralOfcirc(phase, twist, ratio, rmin, ncircles):
r=rmin
g=1+ratio
sc=Graphics()
for _ in range(ncircles):
sc += circle((xc,yc),r,fill=True, facecolor='white')

theta +=twist

xc += g*r*cos(theta); yc +=g*r*sin(theta)
r *= ratio

return sc

show(spiralOfcirc(0,pi/6,1.1,1,12))

r=1.01
g=(1+sqrt(5))/2
ph=2*pi/g
flower=[spiralOfcirc(t*ph,pi/12,r^(t+1),1,12) for t in range(12)]
spirals=Graphics()
for s in flower:
spirals += s
spirals.show(axes=false,xmin=-30,xmax=30,ymin=-30,ymax=30)
save(spirals,"spirals.png",axes=false,xmin=-30,xmax=30,ymin=-30,ymax=30)

sp=spirals(axes=false)

Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python3.9/site-packages/smc_sagews/sage_server.py", line 1234, in execute flags=compile_flags), namespace, locals) File "", line 1, in <module> TypeError: 'Graphics' object is not callable
spirals=Graphics()
n=2
mult=1.1
rmin=.5
ncircs=12
phase=0
turn = 2*pi/12
for k in range(n):

spirals += spiralOfcirc(phase,turn,mult,rmin,ncircs )
phase += 2*turn
#rmin *= mult
spirals.show()







spirs=Graphics()
for f in frames2:
spirs +=f
show(spirs)



nSpirs=3
nTwists=6
ratio=1.05
multiplier=1/ratio

spiralAn=[]
size=5
for tw in range(nTwists):
spiral=Graphics()
#spiral += polygon([(-size,-size), (size,-size), (size,size),(-size,size)])
spiral += circle((0,0),size, fill=False)
for t in range(nSpirs-1):

spiral += spiralOfcirc(tw*2*pi/nTwists+t*2*pi/nSpirs ,pi/12,ratio,r,nSpirs)
r *= multiplier
spiralAn.append(spiral)

frames2=[spiralOfcirc(t*pi/6 ,pi/12,1.2,1,10) for t in range(12)]
spirs=Graphics()
for f in frames2:
spirs +=f
show(spirs)

frames2=[spiralOfcirc(t*2*pi*0.388 ,pi/12,1.01^t,1,10) for t in range(15)]
spirs=Graphics()
for f in frames2:
spirs +=f
show(spirs)

frames2=[spiralOfcirc(t*2*pi*0.388 ,pi/12,1.1,1.05^t,10) for t in range(18)]
spirs=Graphics()
for f in frames2:
spirs +=f
show(spirs)

ph=pi*(1+sqrt(5))
n= 100
mult = 1.00
pic=Graphics()
for t in range(n):
r=sqrt(t)
a = t*ph
#pic += point((r*cos(a),r*sin(a)))
max=20
show(pic,aspect_ratio=1,axes=false,xmin=-max, xmax=max, ymin=-max,ymax=max)
save(pic,"spirals-g.png",aspect_ratio=1,axes=false,xmin=-max, xmax=max, ymin=-max,ymax=max)


show(pic,aspect_ratio=1)



%var r, t
n=13
f=2*pi/13
sp1=Graphics()
for k in range(n):
u=f*k
sp1 += polar_plot(e^(t+u),(t,-u,pi-u),aspect_ratio=1)

n=8
f=2*pi/8
sp2=Graphics()
for k in range(n):
u=f*k
sp2 += polar_plot(e^(pi-t+u),(t,u,pi+u),aspect_ratio=1,color='red')

show(sp2+sp1)



show(sp2+sp1)
save(sp1+sp2,"spir-fib.png",aspect_ratio=1, axes=false)



#transformations needed to build Edmark's tiles

def rot(t,vec):
#rotate vector by angle t counterclockwise
c = cos(t)
s = sin(t)
r= matrix([[c,-s],[s,c]])
return r*vec

def trans(U, vec):
#trans vec by vector U
return vec+U

def cloneDown(A,B,C,D):
#starting with quad ABCD with AB on "bottom", create quad A'B'C'D' below where D'=A, C'=B
h=vector((1,0))
DC = C-D
AB = B-A
r = AB.norm()/DC.norm()  #dilation ratio
U=r*(A-D)
#dilate, centered at A
A1= A
B1= r*(B-A) +A
C1 = r*(C-A)+A
D1=r*(D-A)+A
# now translate so D1 moves to A
A2=trans(U,A1)
B2=trans(U,B1)
C2=trans(U,C1)
D2=trans(U,D1)
#now rotate by t = angle to move CD to be parallel with AB
#we will need a sign of a determinant to know direction
m=Matrix([DC,AB])
t=sgn(m.det())*acos(DC.dot_product(AB)/(DC.norm()*AB.norm()))
A3=rot(t, A2-A)+A
B3=rot(t, B2-A)+A
C3=rot(t, C2-A)+A
D3=rot(t, D2-A)+A

return A3,B3,C3,D3

def cloneLeft(A,B,C,D):
#starting with quad ABCD with AB on "bottom", create quad A'B'C'D' on left where B'=A, C'=D, A'D' on left
X, Y, Z, W = cloneDown(D, A, B, C)
return Y, Z, W, X

def cloneRight(A,B,C,D):
#starting with quad ABCD with AB on "bottom", create quad A'B'C'D' on right where A'=B, D'=C, B'C' on right
X, Y, Z, W = cloneDown(B,C,D,A)
return W,X,Y,Z

#test run, using the coord of the "mother tile"
A=vector((0,0))
B=vector((1,0))
C= vector((1.1056158046757, 1.3611084974398))
D =vector((-0.294552318493650,1.0533128073132))

pic = Graphics()

pic += line([A,B,C,D,A])
X,Y,Z,W = cloneDown(A,B,C,D)
pic += line([X,Y,Z,W,X])
X,Y,Z,W = cloneLeft(A,B,C,D)
pic += line([X,Y,Z,W,X])

X,Y,Z,W = cloneRight(A,B,C,D)
pic += line([X,Y,Z,W,X])

show(pic,aspect_ratio=1)



#draw full tiling

A=vector((0,0))
B=vector((1,0))
C= vector((1.1056158046757, 1.3611084974398))
D =vector((-0.294552318493650,1.0533128073132))

#the four coords below "almost work"
#A=vector((0,0))
#B=vector((1,0))
#C= vector((1, 1.5))
#D =vector((-.5,1))

nsteps =21 #how many levels

pic = Graphics()
pic += line([A,B,C,D,A])
pic += polygon([A,B,C,D],color=hue(1)) #mother quad (so filled in)
for row in range(nsteps):

X,Y,Z,W = cloneDown(X,Y,Z,W)
#pic += line([X,Y,Z,W,X])
pic += polygon([X,Y,Z,W], color=hue(.5*int(mod(row+col+1,2))))
for row in range(nsteps):
for col in range(1, row+3):
X,Y,Z,W = cloneLeft(X,Y,Z,W)
#pic += line([X,Y,Z,W,X])
pic += polygon([X,Y,Z,W], color=hue(.5*int(mod(row+col,2))))

for row in range(nsteps):
for col in range(1, 4):
X,Y,Z,W = cloneRight(X,Y,Z,W)
#pic += line([X,Y,Z,W,X])
pic += polygon([X,Y,Z,W], color=hue(.5*int(mod(row+col,2))))

show(pic,aspect_ratio=1,axes=false)
save(pic,"edmark.pdf",aspect_ratio=1,axes=false)

row

0
.5*int(mod(5,2))

0.500000000000000
a=5
n=12
data =[]
for k in range(n):
print k, a
data.append([k,a])
a = a^2

latex(table(data))

0 5 1 25 2 625 3 390625 4 152587890625 5 23283064365386962890625 6 542101086242752217003726400434970855712890625 7 293873587705571876992184134305561419454666389193021880377187926569604314863681793212890625 8 86361685550944446253863518628003995711160003644362813850237034701685918031624270579715075034722882265605472939461496635969950989468319466936530037770580747746862471103668212890625 9 7458340731200206743290965315462933837376471534600406894271518333206278385070118304936174890400427803361511603255836101453412728095225302660486164829592084691481260792318781377495204074266435262941446554365063914765414217260588507120031686823003222742297563699265350215337206058336516628646003612927433551846968657326499008153319891789578832685947418212890625 10 55626846462680034577255817933310101605480399511558295763833185422180110870347954896357078975312775514101683493275895275128810854038836502721400309634442970528269449838300058261990253686064590901798039126173562593355209381270166265416453973718012279499214790991212515897719252957621869994522193843748736289511290126272884996414561770466127838448395124802899527144151299810833802858809753719892490239782222290074816037776586657834841586939662825734294051183140794537141608771803070715941051121170285190347786926570042246331102750604036185540464179153763503857127117918822547579033069472418242684328083352174724579376695971173152319349449321466491373527284227385153411689217559966957882267024615430273115634918212890625 11 3094346047382578275480183369971197853892556303884969045954098458217021364691229981426352946556259795253241437925401811752196686587974858300172368368138733125186061074284643736990470985187297054554299280845687415532065869107175268273614914079919451498396758201719634752768716320786430383849047583971326289816201205757042613947819180980672005101042630551776963848492580515763228422194205105528219980245495047005616012864622912201600352471713015158045534728307404176895018366960267524059270145304825352506681132363099774878756679292686131556110845011043463378706055597131761908315431198704846311568948881773397779068614881830957489568601480084153293693894869108316525162755529109279622020185143950303962042574196734629855975332530865630162585454286462276973637232439143739157810226810193065511912215461090477312372972236985158496357200031595045778223229181777984057608080727232899920794990550253812318348125460488862923975738912839914768937870513971019200450747608740083773989106089243768942044904738561199032074943292242626646859106725039882814425691037346493790407231123914098041407526301841706149525214943424226821613989004317509345843872795604700727869464468991967905694380830535928512027512496833244483273618867296842941691082751697737211017236389627905793492664437797515223884794021082891431511614939687909024149023453904144326123741074082132970518886413995393818295837875447436245008301682057894055333235883153975009918212890625 \begin{tabular}{ll} $0$ & $5$ \\ $1$ & $25$ \\ $2$ & $625$ \\ $3$ & $390625$ \\ $4$ & $152587890625$ \\ $5$ & $23283064365386962890625$ \\ $6$ & $542101086242752217003726400434970855712890625$ \\ $7$ & $293873587705571876992184134305561419454666389193021880377187926569604314863681793212890625$ \\ $8$ & $86361685550944446253863518628003995711160003644362813850237034701685918031624270579715075034722882265605472939461496635969950989468319466936530037770580747746862471103668212890625$ \\ $9$ & $7458340731200206743290965315462933837376471534600406894271518333206278385070118304936174890400427803361511603255836101453412728095225302660486164829592084691481260792318781377495204074266435262941446554365063914765414217260588507120031686823003222742297563699265350215337206058336516628646003612927433551846968657326499008153319891789578832685947418212890625$ \\ $10$ & $55626846462680034577255817933310101605480399511558295763833185422180110870347954896357078975312775514101683493275895275128810854038836502721400309634442970528269449838300058261990253686064590901798039126173562593355209381270166265416453973718012279499214790991212515897719252957621869994522193843748736289511290126272884996414561770466127838448395124802899527144151299810833802858809753719892490239782222290074816037776586657834841586939662825734294051183140794537141608771803070715941051121170285190347786926570042246331102750604036185540464179153763503857127117918822547579033069472418242684328083352174724579376695971173152319349449321466491373527284227385153411689217559966957882267024615430273115634918212890625$ \\ $11$ & $3094346047382578275480183369971197853892556303884969045954098458217021364691229981426352946556259795253241437925401811752196686587974858300172368368138733125186061074284643736990470985187297054554299280845687415532065869107175268273614914079919451498396758201719634752768716320786430383849047583971326289816201205757042613947819180980672005101042630551776963848492580515763228422194205105528219980245495047005616012864622912201600352471713015158045534728307404176895018366960267524059270145304825352506681132363099774878756679292686131556110845011043463378706055597131761908315431198704846311568948881773397779068614881830957489568601480084153293693894869108316525162755529109279622020185143950303962042574196734629855975332530865630162585454286462276973637232439143739157810226810193065511912215461090477312372972236985158496357200031595045778223229181777984057608080727232899920794990550253812318348125460488862923975738912839914768937870513971019200450747608740083773989106089243768942044904738561199032074943292242626646859106725039882814425691037346493790407231123914098041407526301841706149525214943424226821613989004317509345843872795604700727869464468991967905694380830535928512027512496833244483273618867296842941691082751697737211017236389627905793492664437797515223884794021082891431511614939687909024149023453904144326123741074082132970518886413995393818295837875447436245008301682057894055333235883153975009918212890625$ \\ \end{tabular}
a=5
n=8
data =[]
for k in range(n):
print k, a^2-a
data.append([k,a^2-a])
a = a^2

latex(table(data))

0 20 1 600 2 390000 3 152587500000 4 23283064365234375000000 5 542101086242752217003703117370605468750000000 6 293873587705571876992184134305561419454666388650920794134435709565877914428710937500000000 7 86361685550944446253863518628003995711160003644362813850237034701685918031624270579715074740849294560033595947277362330408531534801930273914649660582654178142547607421875000000000 \begin{tabular}{ll} $0$ & $20$ \\ $1$ & $600$ \\ $2$ & $390000$ \\ $3$ & $152587500000$ \\ $4$ & $23283064365234375000000$ \\ $5$ & $542101086242752217003703117370605468750000000$ \\ $6$ & $293873587705571876992184134305561419454666388650920794134435709565877914428710937500000000$ \\ $7$ & $86361685550944446253863518628003995711160003644362813850237034701685918031624270579715074740849294560033595947277362330408531534801930273914649660582654178142547607421875000000000$ \\ \end{tabular}
# find x,y such that x^2+y^2=p for primes p
def sumOfSquares(n):
# find x,y s.t. x^2+y^=n or if fail, print "no"
output = "no"
u=floor(sqrt(n))
for x in range(u+1):
y= n-x^2
if y.is_square():
output = (x,sqrt(y))
return output
#start value

p=2

#end
biggest = 1000
while p<biggest:
print p, sumOfSquares(p)

#iterate to next prime
p = p.next_prime()


2 (1, 1) 3 no 5 (2, 1) 7 no 11 no 13 (3, 2) 17 (4, 1) 19 no 23 no 29 (5, 2) 31 no 37 (6, 1) 41 (5, 4) 43 no 47 no 53 (7, 2) 59 no 61 (6, 5) 67 no 71 no 73 (8, 3) 79 no 83 no 89 (8, 5) 97 (9, 4) 101 (10, 1) 103 no 107 no 109 (10, 3) 113 (8, 7) 127 no 131 no 137 (11, 4) 139 no 149 (10, 7) 151 no 157 (11, 6) 163 no 167 no 173 (13, 2) 179 no 181 (10, 9) 191 no 193 (12, 7) 197 (14, 1) 199 no 211 no 223 no 227 no 229 (15, 2) 233 (13, 8) 239 no 241 (15, 4) 251 no 257 (16, 1) 263 no 269 (13, 10) 271 no 277 (14, 9) 281 (16, 5) 283 no 293 (17, 2) 307 no 311 no 313 (13, 12) 317 (14, 11) 331 no 337 (16, 9) 347 no 349 (18, 5) 353 (17, 8) 359 no 367 no 373 (18, 7) 379 no 383 no 389 (17, 10) 397 (19, 6) 401 (20, 1) 409 (20, 3) 419 no 421 (15, 14) 431 no 433 (17, 12) 439 no 443 no 449 (20, 7) 457 (21, 4) 461 (19, 10) 463 no 467 no 479 no 487 no 491 no 499 no 503 no 509 (22, 5) 521 (20, 11) 523 no 541 (21, 10) 547 no 557 (19, 14) 563 no 569 (20, 13) 571 no 577 (24, 1) 587 no 593 (23, 8) 599 no 601 (24, 5) 607 no 613 (18, 17) 617 (19, 16) 619 no 631 no 641 (25, 4) 643 no 647 no 653 (22, 13) 659 no 661 (25, 6) 673 (23, 12) 677 (26, 1) 683 no 691 no 701 (26, 5) 709 (22, 15) 719 no 727 no 733 (27, 2) 739 no 743 no 751 no 757 (26, 9) 761 (20, 19) 769 (25, 12) 773 (22, 17) 787 no 797 (26, 11) 809 (28, 5) 811 no 821 (25, 14) 823 no 827 no 829 (27, 10) 839 no 853 (23, 18) 857 (29, 4) 859 no 863 no 877 (29, 6) 881 (25, 16) 883 no 887 no 907 no 911 no 919 no 929 (23, 20) 937 (24, 19) 941 (29, 10) 947 no 953 (28, 13) 967 no 971 no 977 (31, 4) 983 no 991 no 997 (31, 6)

def sumOfSquares3(n):
# find x,y,z such that x^2+y^2+z^2= n  or if fail, print "no"
output = "no"
u=floor(sqrt(n))
for x in range(u+1):
d = n-x^2
if sumOfSquares(d) != "no":
y,z = sumOfSquares(d)
output = x,y,z
return output

for k in [4..100]:
print k, sumOfSquares3(k)

4 (2, 0, 0) 5 (2, 1, 0) 6 (2, 1, 1) 7 no 8 (2, 2, 0) 9 (3, 0, 0) 10 (3, 1, 0) 11 (3, 1, 1) 12 (2, 2, 2) 13 (3, 2, 0) 14 (3, 2, 1) 15 no 16 (4, 0, 0) 17 (4, 1, 0) 18 (4, 1, 1) 19 (3, 3, 1) 20 (4, 2, 0) 21 (4, 2, 1) 22 (3, 3, 2) 23 no 24 (4, 2, 2) 25 (5, 0, 0) 26 (5, 1, 0) 27 (5, 1, 1) 28 no 29 (5, 2, 0) 30 (5, 2, 1) 31 no 32 (4, 4, 0) 33 (5, 2, 2) 34 (5, 3, 0) 35 (5, 3, 1) 36 (6, 0, 0) 37 (6, 1, 0) 38 (6, 1, 1) 39 no 40 (6, 2, 0) 41 (6, 2, 1) 42 (5, 4, 1) 43 (5, 3, 3) 44 (6, 2, 2) 45 (6, 3, 0) 46 (6, 3, 1) 47 no 48 (4, 4, 4) 49 (7, 0, 0) 50 (7, 1, 0) 51 (7, 1, 1) 52 (6, 4, 0) 53 (7, 2, 0) 54 (7, 2, 1) 55 no 56 (6, 4, 2) 57 (7, 2, 2) 58 (7, 3, 0) 59 (7, 3, 1) 60 no 61 (6, 5, 0) 62 (7, 3, 2) 63 no 64 (8, 0, 0) 65 (8, 1, 0) 66 (8, 1, 1) 67 (7, 3, 3) 68 (8, 2, 0) 69 (8, 2, 1) 70 (6, 5, 3) 71 no 72 (8, 2, 2) 73 (8, 3, 0) 74 (8, 3, 1) 75 (7, 5, 1) 76 (6, 6, 2) 77 (8, 3, 2) 78 (7, 5, 2) 79 no 80 (8, 4, 0) 81 (9, 0, 0) 82 (9, 1, 0) 83 (9, 1, 1) 84 (8, 4, 2) 85 (9, 2, 0) 86 (9, 2, 1) 87 no 88 (6, 6, 4) 89 (9, 2, 2) 90 (9, 3, 0) 91 (9, 3, 1) 92 no 93 (8, 5, 2) 94 (9, 3, 2) 95 no 96 (8, 4, 4) 97 (9, 4, 0) 98 (9, 4, 1) 99 (9, 3, 3) 100 (10, 0, 0)
%var x,y
plot3d((x-y)^2+(sqrt(2-x^2)-9/y)^2,(x,0,sqrt(2)),(y,2,3))

3D rendering not yet implemented