/usr/share/genius/examples/explicit-fdm-heat.gel is in genius-common 1.0.23-3.
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 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 | # Category: Differential Equations
# Name: Heat equation solved using Explicit Finite Difference Method (surface plot)
#
# 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
# Boundary conditions by default are
# The boundary conditions are insulated on x=0, u_x(0,t) = 0
# and Dirichlet condition on x=1, u(1,t) = 0
# But this can be changed below.
# Initial value, u(x,0) = ?
function initialf(x) = sin(2*pi*x);
#function initialf(x) = 1.0;
#function initialf(x) = 0.0;
#function initialf(x) = x;
#function initialf(x) = 1-x;
#function initialf(x) = if x < 0.5 then 2*x else 2-2*x
#function initialf(x) = if x < 0.5 then 4*x^2 else 4*(x-1)^2;
#function initialf(x) = if x < 0.5 then -1 else 1;
# Left hand (x=0) endpoint insulated by default
leftendinsulated := true;
leftenddirichlet := 0.0;
# Right hand (x=1) endpoint Dirichlet condition set at 0
# by default
rightendinsulated := false;
rightenddirichlet := 0.0;
# The number of intervals on the x axis
n := 20;
h := float(1/n);
# The value of k/h^2 must be less than or equal 0.5 for the method to be
# stable. We use 0.4 to ensure even the step function is stable.
# Try using 0.55 instead of 0.4, or even try 0.5 with the step function.
# Or set it to less. See what it does to accuracy.
k := 0.4*h^2;
# Maximum t to go up to
maxt := 0.05;
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, set plotevery
# to 1
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))
);
if leftendinsulated then (
u@(1) := u@(2)
) else (
u@(1) := leftenddirichlet
);
if rightendinsulated then (
u@(n+1) := u@(n)
) else (
u@(n+1) := rightenddirichlet
);
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));
|