Make a model in Python using the Forward Euler scheme (FTCS)

Question

Solve the 1-dimensional advection equation ∂θ/dt = -c dθ/dx by using centered differences in space and forward differences in time. Apply a periodic boundary condition θ^(n) (J) = θ^(n) (0). J = 100 and the initial condition is θ^(0)_j = f0(j), where

F0(j)=0, j = 1,...46

F0(47)=0.3

F0(48)=0.6

F0(49)=09

F0(50)= 1.2

F0(51) =0.9

F0(52)=0.6

F0(53)=0.3

F0(54-->)=0

c = 1 #wavenumber

∆x = L/J = 1

∆t = 0.75 #timestep

with this information, make a FTCS scheme model and integrate the model 20 time steps. Plot the solution (numerical and analytical) for the times t=0, and 4 ∆t.

c=1
J=100
L=100
dx=1
dt=0.75

Courant_number = (c*dt)/(2*dx)
F0_j = np.zeros((100,100)) #matrise

F0_j[47,0]=0.3
F0_j[48,0]=0.6
F0_j[49,0]=0.9
F0_j[50,0]=1.2
F0_j[51,0]=0.9
F0_j[52,0]=0.6
F0_j[53,0]=0.3

for j in range(J):
    for n in range(L):
        if n == 99:
            continue
        if j == 99:
            F0_j[j,n] = F0_j[0,n]
            continue
        F0_j[j,n+1] = F0_j[j,n]-(Courant_number*(F0_j[j+1,n]-F0_j[j-1,n]))

plt.plot(F0_j)

Answer 1

On the periodicity condition: With the central difference quotient you propagate the range F0[1:-1,j]. The outer elements need to be copied from the other end. The element 0 gets copied from place J, which is thus not the outer element, and the element J+1 gets copied from place 1. So

F0 = np.zeros((J+2,L+1))  # space, time

Explicit loops are slow, thus should be avoided. Especially as you are already using numpy.

for n in range(L):
    F0[1:-1,n+1] = F0[1:-1,n]-(Courant_number*(F0[2:,n]-F0_j[:-2,n]))
    F0[0,n+1], F0[-1,n+1] = F0[-2,n+1], F0[1,n+1]

You could also use the np.roll procedure to realize this without border cells, C*( np.roll(F0[:,n],-1) - np.roll(F0[:,n],1) )

---

To get a gallery of plots you need to generate each graph separately, either using subplots for a selection of graphs arranged in one figure, or generating multiple figures.

for k in range(11):
    plt.plot(F0[:,10*k], label=f"t={10*k}")
    plt.grid(); plt.legend(); plt.show()