Why does the same algorithm result in different outputs in C++ & Python?

I am running a small code in which there are periodic boundary conditions i.e.,for point 0 the left point is the last point and for the last point zeroth point is the right point. When I run the same code in Python and C++, the answer I am getting is very different.

Python Code

    import numpy as np
    c= [0.467894,0.5134679,0.5123,0.476489,0.499764,0.564578]
    n= len(c)
    Δx= 1.0
    A= 1.0
    M = 1.0
    K=1.0
    Δt= 0.05

    def delf(comp):
        ans = 2*A*comp*(1-comp)*(1-2*comp)
        return ans

    def lap(e,v,w):
        laplace = (w -2*v+ e) / (Δx**2)
        return laplace

    for t in range(1000000):

        μ = np.zeros(n)
        for i in range(n):
            ans1= delf(c[i])
    
            east= i+1
            west= i-1
    
            if i==0:
                west= i-1+n
            if i== (n-1):
               east= i+1-n
        
            ans2= lap(c[east], c[i], c[west])
    
            μ[i] = ans1 - K* ans2

        dc_dt= np.zeros(n)
        for j in range(n):
            east= j+1
            west= j-1
    
            if j==0:
               west= j-1+n
            if j== (n-1):
               east= j+1-n
        
            ans= lap(μ[east], μ[j], μ[west])
    
         dc_dt[j] = M* ans
    
         for p in range(n):
             c[p] = c[p] + Δt * dc_dt[p]

    print(c)

The output in Python 3.7.6 version is

    [0.5057488166795907, 0.5057488166581386, 0.5057488166452102,
     0.5057488166537337, 0.5057488166751858, 0.5057488166881142]

C++ Code

    #include <iostream>

    using namespace std ;

    const float dx =1.0;
    const float dt =0.05;
    const float A= 1.0;
    const float M=1.0;
    const float K=1.0;
    const int n = 6 ;

    float delf(float comp){
        float answer = 0.0;
        answer= 2 * A* comp * (1-comp) * (1-2*comp);
        return answer; }

    float lap(float e, float v, float w){
        float laplacian= 0.0 ;
        laplacian = (e - 2*v +w) / (dx *dx);
        return laplacian; }

    int main(){
        float c[n]= {0.467894,0.5134679,0.5123,0.476489,0.499764,0.564578};
        
        for(int t=0; t<50;++t){
    
            float mu[n];
            for(int k=0; k<n; ++k){
                int east, west =0 ;
                float ans1,ans2 = 0;
        
                ans1= delf(c[k]);
        
                if (k ==0){ west = k-1+n; }
                else{ west = k-1; }
        
                if (k== (n-1)) { east = k+1-n;  }
                else{ east= k+1; }
        
                ans2= lap(c[east], c[k], c[west]);
        
                mu[k] = ans1 - K*ans2; 
            }
        
            float dc_dt[n];
            for(int p=0; p<n; ++p){
                float ans3= 0;
                int east, west =0 ;
        
                if (p ==0){ west = p-1+n; }
                else{ west = p-1;}
        
                if (p== (n-1)) { east = p+1-n; }
                else{ east= p+1; }
        
                ans3= lap(mu[east], mu[p], mu[west]);
        
                dc_dt[p] = M* ans3;
            }

            for(int q=0; q<n; ++q){         
                c[q] = c[q] + dt* dc_dt[q]; 
            }             
        }
        for(int i=0; i<n; ++i){
            cout<< c[i]<< " ";
        }            
        return 0;  
   }

Output in C++ is

    0.506677 0.504968 0.50404 0.50482 0.506528 0.507457

When I am iterating for small steps say t<1000 there is no significant difference in outputs but I am supposed to do this calculation for large number of iterations (in order of 10^7) and here the difference in output is very large.

Answer

I took your code, added the missing closing bracket of the large “for” loop and also changed the length from “50” to “1000000” as in the python version.

Then I replaced all “float” with “double” and the resulting output is:

0.505749 0.505749 0.505749 0.505749 0.505749 0.505749

Thus, of course, implementing the same code in python and in c++ gives the same result. However, the types are obviously important. For example, integers are implemented in a very very different way in python3 than in c++ or almost any other language. But here it is much simpler. Python3 “float” is a “double” in c++ by definition. See https://docs.python.org/3/library/stdtypes.html

Fun fact

the simplest python program that you will have major trouble to reproduce in C++ or most other languages is something like

 myInt=10000000000000000000000000000000345365753523466666
 myInt = myInt*13 + 1
 print (myInt)

since python can work with arbitrary large integers (until your entire computer memory is filled). The corresponding

#include <iostream>
int main(){ 
  long int myInt = 10000000000000000000000000000000345365753523466666;
  myInt = myInt*13 + 1;
  std::cout << myInt << std::endl;
  return 0;
}

will simply report warning: integer constant is too large for its type and will overflow and print a wrong result.