Python

Given this Matlab Code created by my teacher:

function [] = explicitWave(T,L,N,J)
% Explicit method for the wave eq.
% T: Length time-interval
% L: Length x-interval
% N: Number of time-intervals
% J: Number of x-intervals
k=T/N;
h=L/J;
r=(k*k)/(h*h);
k/h

x=linspace(0,L,J+1); % number of points = number of intervals + 1

uOldOld=f(x); % solution two time-steps backwards. Initial condition
disp(uOldOld)
uOld=zeros(1,length(x)); % solution at previuos time-step
uNext=zeros(1,length(x));

% First time-step
for j=2:J
    uOld(j)=(1-r)*f(x(j))+r/2*(f(x(j+1))+f(x(j-1)))+k*g(x(j));
end

% Remaining time-steps
for n=0:N-1
    for j=2:J
        uNext(j)=2*(1-r)*uOld(j)+r*(uOld(j+1)+uOld(j-1))-uOldOld(j);
    end
    uOldOld=uOld;
    uOld=uNext;
end
 plot(x,uNext,'r')
end

I tried to implement this in Python by using this code:

import numpy as np
import matplotlib.pyplot as plt

def explicit_wave(f, g, T, L, N, J):
    """

    :param T: Length of Time Interval
    :param L: Length of X-interval
    :param N: Number of time intervals
    :param J: Number of X-intervals
    :return:
    """
    k = T/N
    h = L/J
    r = (k**2) / (h**2)

    x = np.linspace(0, L, J+1)
    Uoldold = f(x)
    Uold = np.zeros(len(x))
    Unext = np.zeros(len(x))

    for j in range(1, J):
        Uold[j] = (1-r)*f(x[j]) + (r/2)*(f(x[j+1]) + f(x[j-1])) + k*g(x[j])

    for n in range(N-1):
        for j in range(1, J):
            Unext[j] = 2*(1-r) * Uold[j]+r*(Uold[j+1]+Uold[j-1]) - Uoldold[j]

        Uoldold = Uold
        Uold = Unext
    
    
    plt.plot(x, Unext)
    plt.show()

    return Unext, x

However when I run the code with the same inputs, I get different results when plotting them. My inputs:

g = lambda x: -np.sin(2*np.pi*x)
f = lambda x: 2*np.sin(np.pi*x)
T = 8.0
L = 1.0
J = 60
N = 480

Python plot result compared to exact result. The x-es represent the actual solution, and the red line is the function:

Matlab plot result , x-es represent the exact solution and the red line is the function:

Could you see any obvious errors I might have made when translating this code?

In case anyone needs the exact solution:

exact = lambda x,t: 2*np.sin(np.pi*x)*np.cos(np.pi*t) - (1/(2*np.pi))*np.sin(2*np.pi*x)*np.sin(2*np.pi*t)

Answer

I found the error through debugging. The main problem here is the code:

Uoldold = Uold
Uold = Unext

So in Python when you define a new variable as equal to an older variable, they become references to each other (i.e dependent on each other). Let me illustrate this as an example consisting of lists:

a = [1,2,3,4]
b = a
b[1] = 10
print(a)
>> [1, 10, 3, 4]

So the solution here was to use .copy() Resulting in this:

Uoldold = Uold.copy()
Uold = Unext.copy()

Python implementation of Matlab Code – Finite Difference Method

Advertisement

Answer