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.
Advertisement
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.