U-skate "steers itself to horizontal or vertical"

17 views
Skip to first unread message

Robert Munafo

unread,
Jul 25, 2012, 11:48:18 PM7/25/12
to reaction-...@googlegroups.com
On 7/24/12, Tim Hutton <tim.h...@gmail.com> wrote:
> Yes, it's great for drawing the u-skater. After a bit of practice it
> works almost every time.
>
> One thing it highlights is how the grid still interferes with the
> equations. I tried to get the u-skater going at an arbitrary angle and
> it works for a bit but then steers itself to horizontal or vertical.
> At larger scales it lasts longer but the scale has to be increased by
> something like a factor of 8 before the effect stops being obvious. I
> wonder if Penrose grids are less prone to this problem.

You can get the U-skates to "bend" towards the orthogonal directions
by slightly increasing the k parameter (not too far, or it loses its
stability and lengthens indefinitely).

You can get it to bend towards the diagonals by decreasing k slightly
(again, not too far or the U "dies" by turning into a stationary spot)

Better still, you can reduce the tendency to "bend" in any direction
by computing the Laplacian using a Moore neighborhood with a following
weights:

(l+r+u+d)*2/3 + (ul+ur+dl+dr)/6 - (10/3)*center

as opposed to the standard von Neumann neighbourhood, which is:

(l+r+u+d) - 4*center

This also makes the U-skate more stable in coarser simulations, and is
a good idea for all systems involving a Laplacian, i.e. all systems
with diffusion of any kind. One reference is this paper:

Arad, A Yakhot, G Ben-Dor. A Highly Accurate Numerical Solution of a
Biharmonic Equation. Numer. Meth. PDE, 13, pp. 375-397, 1997.
www.bgu.ac.il/~yakhot/homepage/publications/nmpde_4_97.pdf (on page
379)

--
Robert Munafo -- mrob.com
Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 -
mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com

Tim Hutton

unread,
Oct 1, 2012, 10:24:41 AM10/1/12
to reaction-...@googlegroups.com
(replying to an older thread)

Yes, I remember you telling me about this before. After searching,
most people's nine-point stencil has these weights:
1 4 1
4 20 4 / 6
1 4 1
which is of course the same as yours, just expanded out. There are
others mentioned, e.g.
http://en.wikipedia.org/wiki/Discrete_Laplace_operator#Implementation_in_Image_Processing
but perhaps those are for specialized purposes.

Two things confuse me, though:

1) The central cell is weighted at -10/3, not -4. Why is the matrix
above not divided by 5? This would give the same speed of diffusion.
Compare with what Greg Turk says in his 1991 paper: his generalized
Laplacian (with -1 in the center) is multiplied by four to match the
feature size seen on the square grid. Maybe this doesn't matter?

2) Trying different kernels in Ready I'm actually finding that the
following kernel gives the best results:

1/8 7/8 1/8
7/8 -4 7/8 n8_1over8.png
1/8 7/8 1/8

The evidence is attached - the n8_1over8.png wavefront is closest to circular.

Expressed in these terms, the standard 9-point stencil is given by:

1/5 4/5 1/5
4/5 -4 4/5 and is shown in n8_1over5.png
1/5 4/5 1/5

The other images, 1over2, 1over6 and 1over10 are obvious I hope.
n8_1overINF is the standard 5-point stencil:

0 1 0
1 -4 1 n8_1overINF.png (I just wanted the naming
to be consistent)
0 1 0

The pattern I'm testing on is the Oregonator, starting from a single
pixel. (Set up so that Tab advances by 33000 frames, so all the images
are at the same point.) Hopefully there is nothing about this rule
that means I shouldn't be using it to test our Laplacians?

These results confuse me. Why does the standard 9-point stencil not
give me optimal results? I'll go on and test the u-skater 'bending',
which might conceivably give different results.

Tim

--
Tim Hutton - http://www.sq3.org.uk - http://profiles.google.com/tim.hutton/
n8_1over2.png
n8_1over5.png
n8_1over6.png
n8_1over8.png
n8_1over10.png
n8_1overINF.png
oregonator_test.vti

Robert Munafo

unread,
Oct 1, 2012, 5:45:42 PM10/1/12
to reaction-...@googlegroups.com
Thank you Tim for doing these tests!

I apologize for being so assertive and making it sound like "1:4:-20
is the one true best ratio". I found that Arad et al. paper, and maybe
one or two other similar examples of the ratios 1 : 4 : -20, and I
don't really understand the math well enough to prove or disprove it
myself.

On 10/1/12, Tim Hutton <tim.h...@gmail.com> wrote:
> 1) The central cell is weighted at -10/3, not -4. Why is the matrix
> above not divided by 5? This would give the same speed of diffusion.
> Compare with what Greg Turk says in his 1991 paper: his generalized
> Laplacian (with -1 in the center) is multiplied by four to match the
> feature size seen on the square grid. Maybe this doesn't matter?

There's a scaling factor, which I think you already realize. But just
to answer your question directly:

It's -10/3 to ensure the total adds up to 0. The full formula is
"(l+r+u+d)*2/3 + (ul+ur+dl+dr)/6 - (10/3)*center". So the coefficients
are (2+2+2+2)/3 + (1+1+1+1)/6 - 10/3, which adds up to zero.

If you want to use -1 in the center, that's fine as long as you
multiply all the other coefficients by a factor of 3/10 so the sum
ends up still being 0. Then the formula would be "(l+r+u+d)/5 +
(ul+ur+dl+dr)/20 - center".

Since this affects the diffusion speed, we really want whatever
coefficients give a 1.0 scaling factor, so that the value can be
directly used as the Laplacian term in the formula being computed. I
thought that's what the Arad et al. paper was trying to do.

> 2) Trying different kernels in Ready I'm actually finding that the
> following kernel gives the best results:
>
> 1/8 7/8 1/8
> 7/8 -4 7/8 n8_1over8.png
> 1/8 7/8 1/8
>
> The evidence is attached - the n8_1over8.png wavefront is closest to
> circular.

Yes, I agree -- that looks a lot better than the n8_1over5.png one
that I was proposing.

I think that's what we should do, until someone who understands the
math can help us: just experiment with different ratios to see what
produces the most natural-looking results in simulations.

We should also scale the coefficients to make sure the wavefronts are
spreading at the same rate that they would spread if the von Neumann
stencil were being used (i.e. the n8_1overINF.png version). It seems
like you did that, so that pretty much proves that the center should
always be -4, and the others scaled to match.

> The pattern I'm testing on is the Oregonator, starting from a single
> pixel. (Set up so that Tab advances by 33000 frames, so all the images
> are at the same point.) Hopefully there is nothing about this rule
> that means I shouldn't be using it to test our Laplacians?

The choice of the coefficients for the Laplacian is part of the
algorithm for simulating diffusion, so it should be equally valid for
any R-D system, the heat equation, and anything else that uses the
Laplacian. I used Gray-Scott for testing, and I think there are
similar tests for most other R-D systems.

(To do the circular wavefront test in Gray-Scott, set k=F=0.06,
initialize the whole grid to a=1.0, b=0.0, except for a single pixel
which should be a=b=1.0. The paintbrush and pencil tools worked for
me, or I guess you could hack the VTI or VTU text files. You should
get a spot that grows to fill the whole grid, and I predict that with
your n8_1over8 stencil it would grow into a circular shape. If one
pixel doesn't work, set a small square to a=b=1.0.)

> On 26 July 2012 04:48, Robert Munafo <mro...@gmail.com> wrote:
>> [...]
>> Better still, you can reduce the tendency to "bend" in any direction
>> by computing the Laplacian using a Moore neighborhood with a following
>> weights:
>> (l+r+u+d)*2/3 + (ul+ur+dl+dr)/6 - (10/3)*center
>> [...]
>> Arad, A Yakhot, G Ben-Dor. A Highly Accurate Numerical Solution of a
>> Biharmonic Equation. Numer. Meth. PDE, 13, pp. 375-397, 1997.
>> www.bgu.ac.il/~yakhot/homepage/publications/nmpde_4_97.pdf (on page
>> 379)

Tim Hutton

unread,
Oct 2, 2012, 5:22:10 AM10/2/12
to reaction-...@googlegroups.com
On 1 October 2012 22:45, Robert Munafo <mro...@gmail.com> wrote:
> Thank you Tim for doing these tests!
>
> I apologize for being so assertive and making it sound like "1:4:-20
> is the one true best ratio". I found that Arad et al. paper, and maybe
> one or two other similar examples of the ratios 1 : 4 : -20, and I
> don't really understand the math well enough to prove or disprove it
> myself.

No need to apologize - everyone seems to use the 1:4:20 stencil.
It works with one one pixel but doesn't give a circular shape, see
attached. I tried with k=0, F=0.06 and it is more circular but now I'm
thinking that this is highlighting the the situation in my original
Oregonator experiment - we are wrong in expecting these simple growth
patterns to be circular. Maybe?

We need a mathematician. I'll send up the bat signal.

>
>> On 26 July 2012 04:48, Robert Munafo <mro...@gmail.com> wrote:
>>> [...]
>>> Better still, you can reduce the tendency to "bend" in any direction
>>> by computing the Laplacian using a Moore neighborhood with a following
>>> weights:
>>> (l+r+u+d)*2/3 + (ul+ur+dl+dr)/6 - (10/3)*center
>>> [...]
>>> Arad, A Yakhot, G Ben-Dor. A Highly Accurate Numerical Solution of a
>>> Biharmonic Equation. Numer. Meth. PDE, 13, pp. 375-397, 1997.
>>> www.bgu.ac.il/~yakhot/homepage/publications/nmpde_4_97.pdf (on page
>>> 379)
>
> --
> Robert Munafo -- mrob.com
> Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 -
> mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com



gs_Fk0.06.png
grayscott_stencil_test.vti

Tim Hutton

unread,
Oct 2, 2012, 5:46:39 AM10/2/12
to reaction-...@googlegroups.com
Further experiments. I decided to use the heat equation, as it is as
simple as possible. Now we get completely different results: every
image is more or less circular, even with a 5-point stencil, only the
size is changing.

Details:

Each image heat_n8_1over<X>.png used a stencil like this:

1/X 1-(1/X) 1/X
1-(1/X) -4 1-(1/X)
1/X 1-(1/X) 1/X

The image heat_n8_standard9pt.png used the standard stencil from the papers:

1 4 1
4 -20 4 / 6
1 4 1

which is the same as:

1/6 2/3 1/6
2/3 -10/3 2/3
1/6 2/3 1/6

Note that the 1overINF and standard9pt images are near-identical.

Conclusions:
1) Everyone is using the correct stencil. It has those weights for a
good reason.
2) The non-circular patterns seen in Oregonator and Gray-Scott are
actually something to do with the rule!

I wonder what the implications are for Turk's generalized Laplacians.
I'll ask him.
heat_n8_1over2.png
heat_n8_1over5.png
heat_n8_1over8.png
heat_n8_1overINF.png
heat_n8_standard9pt.png
heat_stencil_test.vti

Robert Munafo

unread,
Oct 2, 2012, 5:12:08 PM10/2/12
to reaction-...@googlegroups.com
On 10/2/12, Tim Hutton <tim.h...@gmail.com> wrote:
> On 1 October 2012 22:45, Robert Munafo <mro...@gmail.com> wrote:
>> [...]
>> (To do the circular wavefront test in Gray-Scott, set k=F=0.06,
>> initialize the whole grid to a=1.0, b=0.0, except for a single pixel
>> which should be a=b=1.0. The paintbrush and pencil tools worked for
>> me, or I guess you could hack the VTI or VTU text files. You should
>> get a spot that grows to fill the whole grid, and I predict that with
>> your n8_1over8 stencil it would grow into a circular shape. If one
>> pixel doesn't work, set a small square to a=b=1.0.)
>
> It works with one one pixel but doesn't give a circular shape, see
> attached.

I tried it again, and it still works for me. I don't get a circle
either, I get a sort of square shape with rounded corners. Here is the
VTI file. Open it into Ready and hit the tab key a few times. I have a
Ready "version 0.4" built on May 1st.

I assume my square shape is because the model resolution is too crude.
It's really hard to take a Ready file and change the model resolution
with also affecting the system parmeters, but I am able to get better
circles. I think that increasing the diffusion rates is similar to
improving the spatial resolution, and you also have to increase the
grid size (which destroys the starting pattern) and change the time
step. Of course, it's no longer simulating the same system once you
change the diffusion rates.

> I tried with k=0, F=0.06 and it is more circular but now I'm
> thinking that this is highlighting the the situation in my original
> Oregonator experiment - we are wrong in expecting these simple growth
> patterns to be circular. Maybe?

Considering symmetry, any phenomenon must lack a preference for any
particular direction. If a phenomenon lacks radial symmetry, then we
would still expect the directional features to appear in "random"
directions from one trial to the next (using "random noise" starting
patterns to make them distinct). Since diffusion is involved, any
non-symmetry in the starting pattern (such as the corners of an
initial single pixel) should get smoothed out. If not, then the system
isn't being simulated accurately enough.

The scientific way to test this is to methodically increase the
resolution of delta x and delta t, which implies also increasing the
number of pixels in the grid and increasing the number of time steps.
Here is an example of me doing that with U-skates:

http://mrob.com/sci/talks/sl10-fig1.png

Because our formulas include a Laplacian, which is a 2nd derivative
with respect to space, each doubling of delta-x resolution requires
improving delta-t resolution by a factor of 4 just to prevent
instability. We need another factor of 2 to accomplish the intended
result of more accurate convergence; as a result the doubled
simulation takes 2^2*4*2=32 times as much computation. In my tests I
went in steps of sqrt(2) with each step in delta x, which means
improving delta t by sqrt(8) = 2 sqrt(2), so the tests took sqrt(32)
times as much computation with each improvement.

I found that the rotational speed of a pair of spots drops with each
increase in resolution, and plotting the data showed an asymptotic
tendency towards zero, converging exponentially. By "exponential" I
mean that with each doubling of the number of pixels in the grid, the
time it took for a calibrated change in angle (roughly 20 degrees)
increased by an amount that was roughly constant: so as the number of
pixels increases exponentially, the speed of the rotation drops
exponentially, as measured in the model time scale (not simulation
"time steps", which are unrelated because we'r changing delta t)

I found several other parameters that could be measured (such as the
mass of the U-skate's content of chemical "a" as compared to the
steady-state "background" level, again measured in model units) and
these also had an exponential convergence.

This type of test is currently impossible in Ready, because Ready does
not allow direct manipulation of delta-x and delta-t independently of
the other model parameters (like the diffusion rates of the chemical
system). This makes it hard to compare Ready settings with the results
reported in scientific papers.
grayscott_kf_006_1_pixel.vti

Tim Hutton

unread,
Oct 2, 2012, 6:08:27 PM10/2/12
to reaction-...@googlegroups.com
We can include the delta_x in the formula - attached.
munafo-uskate_with_delta_x.vti

Robert Munafo

unread,
Oct 2, 2012, 7:05:03 PM10/2/12
to reaction-...@googlegroups.com
On 10/2/12, Tim Hutton <tim.h...@gmail.com> wrote:
> We can include the delta_x in the formula - attached.

Oh yes, thank you. I remember now that you did that before. I don't
easily remember how to change the formulas, which is why I made my
program implement all the necessary parameters from the start.

There is still no way to re-size an existing pattern to agree with
changes in delta x. For example, if I take the file you sent, and let
the U-skate run a while so it assumes its characteristic shape, and
then try to edit the dimensions to 512x256x1, the pattern is lost (it
returns to the original plain rectangles). You can try the "Save
Pattern..." command to save the current pattern in a file, but it
doesn't fix the problem: when opening the file, the pattern is there,
but then if you change the dimensions, the pattern gets destroyed
anyway. Using the paintbrush to change the pattern has the same
effect: any work you do gets lost when you resize the grid.

Without a way to re-scale the grid and also preserve the current
pattern, it's very hard to do the kinds of tests I'm proposing.

Robert Munafo

unread,
Oct 2, 2012, 7:13:40 PM10/2/12
to reaction-...@googlegroups.com
On 10/2/12, Tim Hutton <tim.h...@gmail.com> wrote:
> Further experiments. I decided to use the heat equation, as it is as
> simple as possible. Now we get completely different results: every
> image is more or less circular, even with a 5-point stencil, only the
> size is changing.

Oh, darn, well I'm wrong again (-: That's why I stick to programming
most of the time.

I suppose with diffusion only, both chemicals remain in the same
proportion with respect to each other when the pattern spreads. But
because the reaction is non-linear (the rate is a*b^2) an infusion of
chemical b has a greater effect on the reaction rate in a given pixel.
So the impact of diffusion on reaction would be greater in the
orthogonal directions where there is more flow. Or something... I
think it would take a lot of math to show what's going on.

>> We need a mathematician. I'll send up the bat signal.

Yep! This is part of why it was hard for me to publish that paper. The
best I can manage is to throw lots of gigaflops at the problem and
apply known numerical methods. The non-circular effects always went
away when I used larger grids and finer time steps.

Tim Hutton

unread,
Oct 4, 2012, 9:50:53 AM10/4/12
to reaction-...@googlegroups.com
On 3 October 2012 00:05, Robert Munafo <mro...@gmail.com> wrote:
> On 10/2/12, Tim Hutton <tim.h...@gmail.com> wrote:
>> We can include the delta_x in the formula - attached.
>
> Oh yes, thank you. I remember now that you did that before. I don't
> easily remember how to change the formulas, which is why I made my
> program implement all the necessary parameters from the start.
>
> There is still no way to re-size an existing pattern to agree with
> changes in delta x. For example, if I take the file you sent, and let
> the U-skate run a while so it assumes its characteristic shape, and
> then try to edit the dimensions to 512x256x1, the pattern is lost (it
> returns to the original plain rectangles). You can try the "Save
> Pattern..." command to save the current pattern in a file, but it
> doesn't fix the problem: when opening the file, the pattern is there,
> but then if you change the dimensions, the pattern gets destroyed
> anyway. Using the paintbrush to change the pattern has the same
> effect: any work you do gets lost when you resize the grid.
>
> Without a way to re-scale the grid and also preserve the current
> pattern, it's very hard to do the kinds of tests I'm proposing.

Yes, this is a known problem. When we have a way to expand the grid
this should be solved.

For now:

Load that file and resize to 512x256x1 (double). Before running, set
delta_x to 0.25 and timestep to 0.1. The u-skater now works.

Likewise for other sizes.

When you've got a grid that is big enough, you can change delta_x at
any time as long as it is changed slowly - the glider reforms quite
readily at the new size, I find.

It's not ideal but hopefully that will help.

I'm thinking it is a good idea to include delta_x in all our rules. It
makes it clearer for the user, rather than for them have to work out
the maths of how it is absorbed into the diffusion coefficients. Do we
think that is a good idea? It makes Ready a bit more technical. Or
maybe I'll include it in a few rules where it is important and leave
the others, and document the difference somewhere.

Robert Munafo

unread,
Oct 4, 2012, 6:01:50 PM10/4/12
to reaction-...@googlegroups.com
On 10/4/12, Tim Hutton <tim.h...@gmail.com> wrote:
> For now:
>
> Load that file and resize to 512x256x1 (double). Before running, set
> delta_x to 0.25 and timestep to 0.1. The u-skater now works.
>
> Likewise for other sizes.

Thanks Tim! For now I'll have to create all my patterns out of
rectangle lists, by editing the VTI files in Emacs, so that this
technique can be used.

> When you've got a grid that is big enough, you can change delta_x at
> any time as long as it is changed slowly - the glider reforms quite
> readily at the new size, I find.

I do that a lot in 2D. The Uskates are pretty hardy critters (-:
However it doesn't work for large patterns, or when alignment is
critical.

> I'm thinking it is a good idea to include delta_x in all our rules. It
> makes it clearer for the user, rather than for them have to work out
> the maths of how it is absorbed into the diffusion coefficients. Do we
> think that is a good idea? It makes Ready a bit more technical. Or
> maybe I'll include it in a few rules where it is important and leave
> the others, and document the difference somewhere.

Since it is so easy to do with changes to the VTI files, I think we
could address this just by including a few example VTI files that
replicate examples from published papers.

For example, Pearson in 1993 used a "system size" of 2.5 on a 256x256
grid, which makes delta x = 2.5 / 256 = 0.009765625. His D_a and D_b
were 0.00002 and 0.00001 respectively, and timestep was 1.0. Here
(attached) is a VTI file called "Pearson 1993 type Theta" that uses
these settings, with k=0.06 and F=0.04.

I haven't edited the description, we'd want to give the author and
title of the paper and the link to Arxiv (
http://arxiv.org/abs/patt-sol/9304003 ) and we also should explain
that the chemical u is called "a" and v is "b".

I think that would be enough. Anyone who is trying to figure out how
to repeat the results of a scientific paper could download the PDF and
they could figure out the rest pretty easily.
Pearson 1993 type Theta.vti

Tim Hutton

unread,
Oct 5, 2012, 7:26:23 AM10/5/12
to reaction-...@googlegroups.com
Thanks for this file! I've committed it, after adding that info.
Reply all
Reply to author
Forward
0 new messages