H flow on tori

Sha Xin Wei
12/90
415-725-3152

Introduction

See comments at end.
Contrast with with 1/H flow on tori,  also a study of induced flow of profile curve.

We avolve tori of evolution, considering the initial profile curve as a zero set of a function 

	u0: RxR -> R

The  deformation can be recast as an evolution equation for

	u : R^2 x I -> R

where 

	u[p,0] = u0[p]

and

	d/dt u = (-1/2) (u1/x + 
			1/Norm2[Grad[u,{x,y}]] (u11 u2^2 - 2 u12 u1 u2 + u22 u1^2))
			
	[Notation, u1 == d/dx (u) , u2 == d/dy (u) , etc.  ]

Definitions: algebra of evolution equation

Norm2[u_] := u . u

Grad[u,x_List] := D[u@@x,#]& /@ x

Simplified evolution equation

evolu = -u22/2 + (u1*u12*u2 + (u2^2*(u22 - u11))/2)/(u1^2 + u2^2) - u1/(2*x);

Derivation

(-1/2) (u1/x + 1/Norm2[Grad[u,{x,y}]] (u11 u2^2 - 2 u12 u1 u2 + u22 u1^2))

Partial difference operators: (from level_set_flow.ma)

symbolic replacement rules

der2diffrules = {
	u1 -> dD[1,umat],
	u2 -> dD[2,umat],
	u11 -> dD[{1,1},umat],
	u12 -> dD[{1,2},umat],
	u21 -> dD[{2,1},umat],
	u22 -> dD[{2,2},umat],
	y -> ymat,
	x -> xmat
};

dD: pdo on solid toroidal domain, wraparound identifies boundaries:

Clear[dD];

dD[n_Integer,mat_] := RotateLeft[mat] - mat /; (TensorRank[mat]==1);

dD[n_Integer,mat_] := Map[RotateLeft,mat,1] - mat /; (n == 1) && (TensorRank[mat]>1);

dD[2, mat_] := Transpose[Map[RotateLeft, Transpose[mat]]] - mat;

dD[n_Integer,mat_] := RotateLeft[mat] - mat /; (n == 3) && (TensorRank[mat]>2);

dD[{n_Integer},mat_] := dD[n,mat];

dD[nl_List,mat_] := dD[Drop[nl,1],dD[nl[[1]],mat]];

dD[nl_List,mat_,h_] := dD[nl,mat]/h^(Length[nl]);

dD[n_Integer,mat_,h_] := dD[n,mat]/h;

difference form of the evolution

Recall the RHS of the evolution equation:

Applying the der2diffrules , we obtain the evolution equation cast in difference form, with elementwise arithmetic:

Clear[ddt];
ddt[umat_,xmat_] := -dD[1, umat]/(2* xmat) - dD[2, dD[2, umat]]/2 + 
  (dD[1, umat]*dD[2, umat]*dD[2, dD[1, umat]] + 
     (dD[2, umat]^2*(-dD[1, dD[1, umat]] + dD[2, dD[2, umat]]))/2)/
   (dD[1, umat]^2 + dD[2, umat]^2);

evolution scheme & plotting function

plotting

Clear[lplot];
lplot[mat_] := 	ListContourPlot[mat, PlotRange->{-0.1,0.1},
	ContourLevels->3, Axes->{0,0}];
lplot[mat_,origin_] := 	ListContourPlot[mat, PlotRange->{-0.1,0.1},
	ContourLevels->3, Axes->origin];

initial level curve

u0r[x_,y_,r_] := (Abs[x]-1)^2 + y^2 - r^2

u0[x_,y_] := u0r[x,y,.95]

u1[x_,y_] := u0r[x,y,.9]

Clear[MakeData];

MakeData[fun_,xrange_,yrange_] := Block[{},
	umat0 = Table[fun, xrange,yrange];
	xmat = Transpose[Table[xrange[[1]],xrange,yrange]];
	ymat = Table[xrange[[1]],xrange,yrange];
	xy0 = {1 - xrange[[2]]/xrange[[4]], 1 - yrange[[2]]/yrange[[4]]};
	Print["Dimensions[umat0] = ", Dimensions[umat0],
	"; Dimensions[xmat] = ",Dimensions[xmat], "; xy0 = ",xy0]
]

evolution

Clear[EvolveCurve];
EvolveCurve[um0_,xs_,dt_,n_] := Block[{um = um0,history},
	history = {};
	Do[Block[{},
		um += ddt[um,xs] dt;
		history = Append[history,um];
		lplot[um]],
	{ki,1, n}
	];
	history
];
EvolveCurve[um0_,xs_,dt_,n_,origin_] := Block[{um = um0,history},
	history = {};
	Do[Block[{},
		um += ddt[um,xs] dt;
		history = Append[history,um];
		lplot[um,origin]],
	{ki,1, n}
	];
	history
]

Experiments

Needs["Utilities`ShowTime`"]

Half of torus:

MakeData[u0[xi,yi], {xi,-.5,2.5,3./39},{yi,-1.5,1.5,3./39}]

xy1 = xy0;
umat1t = Transpose[umat0];
xmat1 = Chop[Abs[xmat]];

lplot[umat1t,xy1]

The lines mark u = -0.1, 0.0, and 0.1; u < 0 in the interior.

uhistory1 = EvolveCurve[umat1t,xmat1, -0.1,25,xy1];

To animate: select the group below by clicking between the righthand vertical lines, and then choose Graph/Animate Selected Graphics from the menu.

525.483 Second on Lasso 39x39 x 25

Comments:

My initial evolution equation was incorrect because I allowed some x<0 when in fact the horizontal -- radial -- coordinate should always be nonnegative.

It seems that a fat (1, 0.9) torus does evolve into a (non-round) sphere, but a torus with radii (1, 0.6) also tends to fill in its central hole.  I'll do a bit more next week.

Another problem is to check that the axial symmetry is indeed preserved as expected.  I'll run another simulation with axially symetric initial data (double circle).