Python implementation of Matlab Code - Finite Difference Method

Question

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)

It is possible to spot the issue by working methodically: debug both implementations and step through them line by line, inspecting the intermediate results at each stage. When they diverge - that's your bug. — Dev-iL
– Dev-iL, Commented Nov 8, 2020 at 13:21
Thanks for your advice, I found the error. I will make an Edit showcasing the error in case someone makes the same mistake as me — Bob Pen
– Bob Pen, Commented Nov 8, 2020 at 14:06

Bob Pen · Accepted Answer · 2020-11-08 14:12:38Z

1

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()

answered Nov 8, 2020 at 14:12

Bob Pen

676 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Dev-iL Over a year ago

Note that this could also happen in MATLAB, when dealing with handle classes\objects.

Collectives™ on Stack Overflow

Python implementation of Matlab Code - Finite Difference Method

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related