Different output for equivalent code between Matlab and Julia

Question

I have a working Matlab algorithm for the evolution of a pulse via solving the 1D wave equation. I know the code works and I can see the pulse traveling to the sides and bouncing off the walls, when plotting; but when I translate the corresponding steps as Julia code, however, the results are completely different for the same set of inputs. I insert a simple pulse that should propagate across a string in both directions but that isn't exactly the result that I'm getting in Julia... and I don't know why.

Could anybody help point out what I'm doing wrong?

Thank you,

This is the corresponding Julia code:

using Plots,OhMyREPL

    #-solving a PDE.
    
    #NOTE: DOMAIN
    # space
    Lx = 20
    dx = 0.1
    nx = Int(trunc(Lx/dx))
    x  = range(0,Lx,nx)#0:dx:Lx-1
    
    
    #NOTE: TIME
    T = 10
    
    #NOTE: Field variable.
    #-variables.
    wn     = zeros(nx)#-current value.
    wnm1   = wn#zeros(nx)#-w at time n-1.
    wnp1   = wn#zeros(nx)#-w at time n+1.
    # wnp1[1:end] = wn[1:end]#-w at time n+1.
    
    #NOTE: PARAMETERS.
    CFL  = 1#-CFL = c.dt/dx
    c   = 1#-the propagation speed.
    dt  = CFL * dx/c    
    
    #     #NOTE: Initial Conditions.
    wn[49] = 0.1; wn[50] = 0.2;  wn[51] = 0.1
    wnp1[1:end] = wn[1:end]#--assumption that they are being equaled.
    wnp1 = wn#--assumption that they are being equaled.
    
    
    run = true
    if run == true
        #NOTE: Initial Conditions.
        wn[49] = 0.1; wn[50] = 0.2;  wn[51] = 0.1
        wnp1 = wn#--wnp1[1:end] = wn[1:end]#--assumption that they are being equaled.
        
        # NOTE: Time stepping loop.
        t = 0
        while t<T
            #-apply reflecting boundary conditions.
            wn[1]=0.0;wn[end]=0.0
    
            #NOTE: solution::
            global t = t+dt
            global wnm1 = wn; global wn = wnp1#-save CURRENT and PREVIOUS arrays.
            for i=2:nx-1
                global wnp1[i] = 2*wn[i] - wnm1[i] + CFL^2 * ( wn[i+1] -2*wn[i] + wn[i-1]  )
            end
    
        end
        plot(x,wnp1, yrange=[-0.2,0.2])
    end

Only one peak

This is the corresponding Matlab code:

%Domain
%Space.
Lx = 20;
dx = 0.1;
nx = fix(Lx/dx);
x  = linspace(0,Lx,nx);

%Time
T = 10;

%Field Variables

wn = zeros(nx,1);
wnm1 = wn;
wnp1 = wn;
%Parameters...
CFL = 1;
c = 1;
dt = CFL*dx/c;


%Initial conditions
wn(49:51) = [0.1 0.2 0.1];    
wnp1(:) = wn(:);




run = true;
if run == true
    %Initial conditions
    wn(49:51) = [0.1 0.2 0.1];    
    wnp1(:) = wn(:);
    
    %Time stepping loop
    t = 0;
    while t<T      
      %Reflecting boundary conditions:
      wn([1 end]) = 0.0;
      %Solution:
      t=t+dt;
      wnm1=wn;wn = wnp1;%-save current and previous arrays.
      for i=2:nx-1
        wnp1(i) = 2*wn(i) - wnm1(i) + CFL^2 * (  wn(i+1) - 2*wn(i) + wn(i-1)   );
      endfor    
    end
    plot(x,wnp1)
end

Two peaks, as it should be.

If anybody could point me in the right direction, I'd greatly appreciate it.

Answer 1

Julia arrays are designed to be fast. So, when an array is assigned to another array, the default behaviour is make the two arrays the same - same variable address, occupying the same area of memory. Editing one array does the exact same edit on the other array:

 x = [1,2,3]
 y = x

 x[1] = 7

 # x is [7,2,3]
 # y is [7,2,3] also

To copy data from one array to another an explicit copy is required:

x = [1,2,3]
y = copy(x)

x[1] = 7

# x is [7,2,3]
# y is [1,2,3] still

The same thing happens with Python lists:

x = [1,2,3]
y = x

x[0] = 7 # Julia indexes arrays 1+, Python indexes lists 0+

# x is [7,2,3]
# y is [7,2,3] also

It is done this way to minimise the effort of unnecessary copying, which can be expensive.

The upshot is that if the code is running but gives the wrong result, a possible cause is assigning an array when a copy is intended.

The question author indicated that this was indeed the cause of the code's behaviour.