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.

CodePudding user response：

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