[genius] Tue Apr 21 00:35:54 2015 Jiri (George) Lebl <jirka 5z com>

[George Lebl <jirka src gnome org>]

commit f5eb07ee34e5cf344fe6b67520848feade8f6417
Author: Jiri (George) Lebl <jiri lebl gmail com>
Date:   Tue Apr 21 00:35:58 2015 -0500

    Tue Apr 21 00:35:54 2015  Jiri (George) Lebl <jirka 5z com>
    
        * examples/laplace-fdm.gel: update to add more things to play around
          with such as nonhomogenity and make it easy to change boundary
          conditions.  Also rename points to intervals as that's what it
          actually is, and document the thing a bit more.

 ChangeLog                |    7 ++++
 examples/laplace-fdm.gel |   69 ++++++++++++++++++++++++++++++++++------------
 2 files changed, 58 insertions(+), 18 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 8b2118a..2819bf6 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+Tue Apr 21 00:35:54 2015  Jiri (George) Lebl <jirka 5z com>
+
+       * examples/laplace-fdm.gel: update to add more things to play around
+         with such as nonhomogenity and make it easy to change boundary
+         conditions.  Also rename points to intervals as that's what it
+         actually is, and document the thing a bit more.
+
 Wed Apr 08 11:39:22 2015  Jiri (George) Lebl <jirka 5z com>
 
        * examples/vibrating-drumhead-modes.gel: make n consistently come
diff --git a/examples/laplace-fdm.gel b/examples/laplace-fdm.gel
index 17efa87..f795d0e 100644
--- a/examples/laplace-fdm.gel
+++ b/examples/laplace-fdm.gel
@@ -1,26 +1,50 @@
 # Category: Differential Equations
 # Name: Laplace equation solution using Finite Difference Method
+#
+# Solve Laplace equation using FDM with Liebmann's iterative method
+# to solve the resulting linear equations.  The animation is really
+# about the iterations of Liebmann's method and how it converges.
+#
+# The equation is Lu = f(x,y) on the square [0,pi]^2.  Normally f(x,y)
+# is zero but you can change it.
 
-points = 25;
+
+# The number of intervals on each side.  There will be (intervals+1)^2
+# total points to solve for and hence that many equations.
+intervals = 20;
 
 # zero initial guess
-u = zeros(points+1,points+1);
+u = zeros(intervals+1,intervals+1);
 
-# how about random initial guess
-#u = 2*rand(points+1,points+1)-ones(points+1,points+1);
+# how about random initial guess (interval [-1,1])
+#u = 2*rand(intervals+1,intervals+1)-ones(intervals+1,intervals+1);
 
-# initial guess of -1s, that's actually somewhat slow to converge
-#u = -ones(points+1,points+1);
+# initial guess of -1s, that's actually somewhat slow to converge if f=0
+#u = -ones(intervals+1,intervals+1);
 
 # The side conditions
-for n=1 to points+1 do
-  u@(n,1) = -sin(pi*(n-1)/points);
-for n=1 to points+1 do
-  u@(n,points+1) = sin(2*pi*(n-1)/points);
-for n=1 to points+1 do
-  u@(1,n) = 0.5*sin(pi*(n-1)/points);
-for n=1 to points+1 do
-  u@(points+1,n) = 0;
+function bottomside(x) = -sin(x);
+function topside(x) = sin(2*x);
+function leftside(y) = 0.5*sin(y);
+function rightside(y) = 0;
+
+# nonhomogeneity, the f in Lu=f
+function f(x,y) = 0;
+#function f(x,y) = 2;
+#function f(x,y) = -2
+#function f(x,y) = 6*exp(-(x-pi/2)^2-(y-pi/2)^2);
+#function f(x,y) = 10*sin((x-pi/2)*(y-pi/2));
+
+
+# set up the side conditions
+for n=1 to intervals+1 do
+  u@(n,1) = bottomside(pi*(n-1)/intervals);
+for n=1 to intervals+1 do
+  u@(n,intervals+1) = topside(pi*(n-1)/intervals);
+for n=1 to intervals+1 do
+  u@(1,n) = leftside(pi*(n-1)/intervals);
+for n=1 to intervals+1 do
+  u@(intervals+1,n) = rightside(pi*(n-1)/intervals);
 
 # don't draw the legend
 SurfacePlotDrawLegends = false;
@@ -29,18 +53,25 @@ PlotWindowPresent(); # Make sure the window is raised
 
 # plot the data
 SurfacePlotDataGrid(u,[0,pi,0,pi]);
+
 # If you want to export the animation to a sequence of .png
 # images uncomment here and below
 #ExportPlot ("animation" + 0 + ".png");
 
 for n = 1 to 300 do (
+  # to slow down the animation
   wait (0.1);
 
   maxch = 0;
-  for i=2 to points do (
-    for j=2 to points do (
-      old = u@(i,j); 
-      u@(i,j) = (1/4)*(u@(i-1,j)+u@(i+1,j)+u@(i,j-1)+u@(i,j+1));
+  for i=2 to intervals do (
+    for j=2 to intervals do (
+      old = u@(i,j);
+      
+      # This is the equation coming from FDM
+      u@(i,j) = (1/4)*(u@(i-1,j)+u@(i+1,j)+u@(i,j-1)+u@(i,j+1))
+        - (1/4)*f(pi*(i-1)/intervals,pi*(j-1)/intervals)*(pi/intervals^2);
+
+      # Keep track of maximim change
       if |u@(i,j)-old| >  maxch then
         maxch = |u@(i,j)-old|
     )
@@ -49,5 +80,7 @@ for n = 1 to 300 do (
 
   # plot the data
   SurfacePlotDataGrid(u,[0,pi,0,pi]);
+
+  # Animation saving
   #ExportPlot ("animation" + n + ".png");
 );