Useful GEL code for teaching (or learning) - Differential Equations

Cut and paste the code into a new program and hit "Run" in the GNOME version.

D'Alembert solution animations

First, an initial pulse on a finite length string with fixed ends

# A periodic pulse function

function F(x) = (

while x < -1 do x = x + 2;

while x > 1 do x = x - 2;

if x < 0 then

-F(-x)

else if x < 0.45 or x > 0.55 then

0

else (

if x < 1/2 then

20*(x-0.45)

else

20*(0.55-x)

)

)

LinePlotDrawLegends = false;

# Run the animation

for n=1 to 1000 do (

t = n*0.005;

LinePlot(`(x)= (F(x-t) + F(x+t))/2,[0,1,-1.1,1.1])

)

Second, an initial square wave with zero initial velocity on an unbounded string (or at least not bounded within [-5,5])

function F(x) = (

if x < -1 or x > 1 then

0

else

1

)

LinePlotDrawLegends = false;

for n=1 to 300 do (

t = n*0.01;

LinePlot(`(x)=(F(x-t)+F(x+t))/2,[-5,5,-0.1,1.1])

)

And lastly, in initial position of 0 but the initial velocity of the square wave on an unbounded string

function G(x) = (

if x < -1 or x > 1 then

0

else

1

)

function integralofG(a,b) = (

if b < a then

integralofG(b,a)

else if b==a or b < -1 or a > 1 then

0

else (

min(b,1) - max(a,-1)

)

)

LinePlotDrawLegends = false;

for n=1 to 600 do (

t = n*0.01;

LinePlot(`(x)=(integralofG(x-t,x+t))/2,[-5,5,-0.1,1.1])

)

Standing waves animation

LinePlotDrawLegends = false;

for n=1 to 1000 do (

t = n*0.015;

LinePlot(`(x)=cos(t)*sin(x),[-5,5,-2,2])

)

Fourier series

Computing the Fourier series for a half sawtooth for 1 to 20 harmonics.

# The Fourier series up to nth harmonic

function ffsn(x,n) = (

(1/4) +

sum k=1 to n do (

(1/(k*pi)^2)*((-1)^k - 1) * cos(k*pi*x) + ((-1)^(k+1)/(k*pi)) * sin(k*pi*x)

)

)

LinePlotWindow=[-2,2,-0.5,1.5];

LinePlotDrawLegends=false;

LinePlotClear();

# Draw the graph by hand to correctly draw the gap

LinePlotDrawLine([-2,0;-1,1],"thickness",2,"color","darkblue");

LinePlotDrawLine([-1,0;0,0;1,1],"thickness",2,"color","darkblue");

LinePlotDrawLine([1,0;2,0],"thickness",2,"color","darkblue");

# Wait for 1 second before starting with the Fourier series

wait(1);

# Do first 20 harmonics

for n=1 to 20 do (

# compute values and put in a vector instead of using the built in function grapher

v=null;

for x=-2 to 2 by 0.003 do

v=[v,[x,ffsn(x,n)]];

LinePlotClear ();

# Draw the graph by hand to correctly draw the gap

LinePlotDrawLine([-2,0;-1,1],"thickness",2,"color","darkblue");

LinePlotDrawLine([-1,0;0,0;1,1],"thickness",2,"color","darkblue");

LinePlotDrawLine([1,0;2,0],"thickness",2,"color","darkblue");

# Plot the fourier series up to nth harmonic

LinePlotDrawLine(v,"thickness",2,"color","darkgreen");

)

Shock waves, traffic flow equation (nonlinear first order PDE)

If you have the equation u_t + g(u) u_x = 0 with initial condition u(x,0)=phi(x), you can solve by finding characteristics, but you can come up with shocks (discontinuities). So here's a way to draw the characteristics to see where shocks form:

# We are looking at characteristics (looking for shocks)

# for the equation u_t + g(u) u_x = 0, with initial condition

# u(x,0) = phi(x)

# this is for flux u*(2-u)

#function g(u) = 2-2*u;

# Quadratic flow for flux u^2 / 2

function g(u) = u;

# ramp down

function phi(x) = (

if x < 0 then 1

else if x < 1 then 1-x

else 0

);

# rampup

#function phi(x) = (

# if x < 0 then 0

# else if x < 1 then x

# else 1

#);

# we'll have x going from

xstart = -3;

xend = 3;

xstep = 0.1;

# end at this t, start at 0

tend = 1.5;

LinePlotWindow = [xstart,xend,-0.1,tend+0.1];

LinePlotDrawLegends = false;

LinePlotClear();

# Draw characteristics

for x0=xstart to xend by xstep do (

LinePlotDrawLine([x0,0,x0+g(phi(x0))*tend,tend])

)

Using RungeKuttaFull to draw solutions of differential equations

If we have just one first order ODE:

# clear the plot