InterdiffusionICS2025/resol1.py at main · ClementPerrochaud/InterdiffusionICS2025 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erf

def RK4(U0,X,f,para):
    lU = np.zeros((len(X),2))
    lU[0] = U0
    print(U0)
    dx = X[1]-X[0]
    for i,x in enumerate(X[:-1]):
        print(lU[i])
        k1 = f(lU[i],         x,      para)
        k2 = f(lU[i]+k1*dx/2, x+dx/2, para)
        k3 = f(lU[i]+k2*dx/2, x+dx/2, para)
        k4 = f(lU[i]+k3*dx,   x+dx,   para)
        lU[i+1] = lU[i] + (k1 + 2*k2 + 2*k3 + k4)*dx/6
    return lU

def f(U,x,para):
    u,v = U
    d,n = para
    return np.array([ -2/d*x*u*v**(-n/(n+1)), u ])

def get_c(d,n,c0,X):
    U0 = np.array([-1e-12,c0**(n+1)])
    para = (d,n)

    # integration
    print("bip")
    U = RK4(U0,X,f,para)
    print("boop")
    v = U[:,1]

    # normalisation
    vmax = np.max(v)
    vmin = np.min(v)
    v = (v-vmin)/(vmax-vmin)

    c = v**(1/(n+1))
    return c

x0 = -5
x1 = 5
nx = 1000

X  = np.linspace(x0,x1,nx)

c_0 = get_c(1,0,1,X)
c_1 = get_c(1,1,1,X)
c_2 = get_c(1,2,1,X)
c_3 = get_c(1,3,1,X)
c_4 = get_c(1,4,1,X)

cref0 = (1+erf(-X))/2

# plot
fig,ax = plt.subplots()
ax.plot(X, cref0, lw=8, color="cyan", alpha=0.6, zorder=9,  label="analytical solution for D=d")
ax.plot(X, c_0,   lw=3, color="darkblue",        zorder=10, label="numerical integration for D=d")
ax.plot(X, c_1,   lw=3, color="darkgreen",       zorder=11, label="numerical integration for D=dc")
ax.plot(X, c_2,   lw=3, color="darkorange",      zorder=12, label="numerical integration for D=dc2")
ax.plot(X, c_3,   lw=3, color="darkred",         zorder=13, label="numerical integration for D=dc3")
ax.plot(X, c_4,   lw=3, color="indigo",      zorder=14, label="numerical integration for D=dc4")
ax.legend()
ax.grid()

plt.show()