/usr/share/genius/examples/explicit-fdm-heat.gel is in genius-common 1.0.21-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 | # Category: Differential Equations
# Name: Heat equation solved using Explicit Finite Difference Method
#
# Solve heat equation using explicit FDM. We try to animate as different
# times are solved for, though of course we must skip graphing some times when
# there are too many steps.
#
# NOTE that this program changes your variable names in the surface plotting.
# To set them back to default do SurfacePlotVariableNames = ["x","y","z"];
#
# The equation is u_t = u_{xx} where 0 < x < 1
# The boundary conditions are insulated on x=0, u_x(0,t) = 0
# and Dirichlet condition on x=1, u(1,t) = 0
# The number of intervals on the x axis
n := 20;
h := float(1/n);
function initialf(x) = sin(2*pi*x);
# The value of k/h^2 must be less than or equal 0.5 for the method to be
# stable. try putting 0.55 in the line below instead of 0.5. Try setting it
# to less than 0.5 to what it does to accuracy
k := 0.5*h^2;
# Maximum t to go up to
maxt := 0.06;
m := round(maxt/k);
# Never plot more than 4*n steps in the t (that is y) direction
# rather skip rows. If you do not want to skip rows, comment out
# the next line.
plotevery := max(1,floor(m/(2*n)));
#Set up initial value
u := null;
data := null;
for j=1 to n+1 do (
x := (j-1)*h;
u@(j) := initialf(x);
data := [data;[x,0,u@(j)]]
);
SurfacePlotDrawLegends := true;
SurfacePlotVariableNames := ["x","t","u"];
PlotWindowPresent(); # Make sure the window is raised
# If you change it from 0 to plotevery, then no plots will be generated
toplot := 0;
#toplot := plotevery;
for i=1 to m do (
v := u;
for j=2 to n do (
u@(j) := v@(j) + (k/(h^2))*(v@(j+1)+v@(j-1)-2*v@(j))
);
u@(1) := u@(2);
u@(n+1) := 0;
increment toplot;
if toplot == plotevery then (
toplot := 0;
for j=2 to n do (
data := [data;[(j-1)*h,i*k,u@(j)]]
);
data := [data;[0,i*k,u@(1)]];
data := [data;[1,i*k,u@(n+1)]];
SurfacePlotData(data,"u",[0,1,0,maxt])
)
);
#print left hand endpoint
print ("Approximate value of u(0," + maxt + ") ~~ " + u@(1));
|