Advanced Search

Question: Iterative method in 3 dimensions

Question:Iterative method in 3 dimensions

Maplehelpee123 55 Maple
iterator
April 20 2013
0

Hi Everyone:)

 

I'm doing a mini project at the moment and I'm pretty stuck with how to go about something on Maple...

 

My aim is to get the catenoid from an initial cylinder by use of an iterative method.

To find the minimal surface r=g(θ,z) - the catenoid - I want to minimise the Dirichlet integral (in cylindrical polar coordinates)

so that for ,

 

and

 

I have the initial points on the cylinder (mostly thanks to people on this website)

 

Numerical Approximation of Minimal Surfaces:

 

 

restart; with(plots); with(plottools); plots[interactive]; with(ArrayTools)

NULL

Getting the points of the cylinder to be moved:

 

 

slices := [seq(polarplot([2, theta, theta = 0 .. 2*Pi], numpoints = 50), z = 0 .. 2, .1)]

slices[10]

op(slices[10]); CURVES(Vector(4, {1 = ` 1..143 x 1..2 `*Array, 2 = `Data Type: `*float[8], 3 = `Storage: `*rectangular, 4 = `Order: `*C_order}), Vector(4, {1 = ` 1..33 x 1..2 `*Array, 2 = `Data Type: `*float[8], 3 = `Storage: `*rectangular, 4 = `Order: `*C_order}), Vector(4, {1 = ` 1..33 x 1..2 `*Array, 2 = `Data Type: `*float[8], 3 = `Storage: `*rectangular, 4 = `Order: `*C_order}), COLOUR(RGB, 1., 0., 0.)), AXESLABELS(x, y)

a := op([1, 1], slices[10])

pointplot(convert(a, listlist))

data := Array(1 .. 1, 1 .. 3); Matrix(1, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0})

for j to nops(slices) do if nops(op(1, slices[j])) > 1 then data := Concatenate(1, data, seq(Concatenate(2, op([1, k], slices[j]), Vector[column]([seq(-.1+.1*j, s = 1 .. op(2, Dimensions(op([1, k], slices[j]))[1]))])), k = 1 .. nops(op(1, slices[j]))-1)) end if end do; data; Vector(4, {1 = ` 1342 x 3 `*Matrix, 2 = `Data Type: `*anything, 3 = `Storage: `*rectangular, 4 = `Order: `*Fortran_order})

points := pointplot3d(data, axes = framed, colour = black)NULL

points

getdata(points)

ThePoints := convert(data, list, nested = true)

T2 := [seq(ThePoints[i], i = 2 .. 51)]NULL

BC1 := pointplot3d(T2, axes = framed, colour = black)

BC1

T3 := [seq(ThePoints[i], i = 1002 .. 1051)]

BC2 := pointplot3d(T3, axes = framed, colour = black)

BC2

cylinderpic := plot3d(2, theta = 0 .. 2*Pi, z = 0 .. 2, coords = cylindrical, axes = framed, color = white, transparency = .75); cylinderpic

T := [seq(ThePoints[i], i = 52 .. 1001)]``

Nodes := pointplot3d(T, axes = framed, colour = cyan)

Nodes

display(BC1, BC2, cylinderpic, Nodes)

  

``

Convert the points in to cylindrical coordinates:

 

 

SortedPts := sort(ThePoints[2 .. -1], proc (A, B) options operator, arrow; evalb(A[1] < B[1] or A[1] = B[1] and A[2] < B[2] or A[1 .. 2] = B[1 .. 2] and A[3] < B[3]) end proc)

CylPts := `~`[proc (P) options operator, arrow; [sqrt(P[1]^2+P[2]^2), arctan(P[2], P[1]), P[3]] end proc](SortedPts)

remove(proc (P) options operator, arrow; fnormal(P[1]-2.0) = 0. end proc, CylPts)

The Matrices R[i,j], θ[i,j] and Z[i,j] for using in the iteration: NULL

 

 

z := Matrix(21, 50, proc (i, j) options operator, arrow; CylPts[21*j-21+i][3] end proc, datatype = float[8]); theta := Matrix(21, 50, proc (i, j) options operator, arrow; CylPts[21*j-21+i][2] end proc, datatype = float[8]); r := Matrix(21, 50, fill = 2., datatype = float[8])

NULL

NULL

 

NULL

 



Download The_Whole_Program_.mw

 

I have the relaxation methods 

a[i]^((k+1))=a[i]^((k))-w(G[a[i]](a[i,k]))/(G[a[i]a[i]](a[i,k]))``

and

Download relaxation_method.mw

 

where a[i] are all of the coordinates of the surface, λ[i] is the maximum eigenvalue of the Hessian and w is the relaxation parameter.

When I differentiated G with respect to a[i], it ended up with 0 in the RHS of the first method already due to the radius values all being the same (cylinder).  The paper I've been looking at agrees with this being zero, so they suggest to use the second relaxation method... Calculating the hessian was fine but the eigenvalues took so long to get out I stopped the calculation because it seemed unlikely to need that long.

If anyone has any ideas on how to do this I'll be eternally grateful!

 

Thanks,

Rach
Please Wait...