冰川长度的数值模拟。编码错误

Question

我目前正在做我的地球物理学学士学位项目，重点是格陵兰冰川的演化。我的项目部分由编程组成，因为我想对此做得更好。但是，我撞墙了。我似乎无法弄清楚我的代码出了什么问题，而且我已经尝试了好几天了。我希望任何输入/反馈。

基本问题是描述冰川如何随时间演化，这是通过使用Eulers方法求解微分方程（dLdt）来完成的。

我已经尝试制作一个简单的for循环，但对于冰川（L）的长度却得到了负值，这是不可能发生的，因此很显然是错误的。我只是不知道在哪里。床配置文件以及所有常量和公式都是正确的，因此我只是以错误的方式循环。

import numpy as np
import matplotlib.pyplot as plt
from math import exp # exponential function
from math import *

n = 5000    #Iterations (5000 years)
t = np.ones(n)
t[0] = 0    #Time(0)
dt = 1      #Time step

    ## Parameters ##
d0 = 200     #d0 (Start height)
s = 0.014   #Slope
lambd = 300     #lambda
sigma = 10000   
xs = 40000  #x-coordinate of the "bump"
ac = 0.0005 #Accumulation
alfaf = 0.7 #Constant
alfam = 2 #Constant
eps = 1     #Constant
delta = 1.127 #Constant
c = 2.4 #Climate sensitivity

    ## Grids ##
x = np.linspace(0,50000,5000) #My x grid
x[0] = 200              #Start position of glacier's head
L = np.ones(n)          #Glacier Length
B = np.zeros(n)
a = np.ones(n)*ac
Hf = np.ones(n)         #Defining vectors to loop through
F = np.ones(n)
dLdt = np.zeros(n)

for i in range(1, n):
    t[i] = i
    a[i] = a[i-1]+1*ac  #Accumulation (constant)
    x[i] = d0-s*x[i]+lambd*np.exp(-((x[i]-xs)/(sigma))**2) #Bed geometry
    B[i] = a[i-1]*L[i-1]    #Surface mass balance
    Hf[i] = np.max([alfaf*np.sqrt(L[i-1]), -eps*delta*x[i-1]]) #Ice thickness at glacier front
    F[i] = np.min([0, c*x[i-1]*Hf[i-1]])    #Flux of ice at the glacier front
    dLdt[i] = (2*(B[i-1]+F[i-1]))/(3*alfam*np.sqrt(L[i-1])) #Differential equation
    L[i] = L[i-1] + 1*dLdt[i-1]     #Updating L with Eulers Method

plt.plot(t,L)

Answer 1

如果最后一个L[i - 1]值小于零，则可以简单地中断循环。不要忘记从数据数组L（使用slicing）中进行绘制（以及从时间数组t将其维数减小为L数组的大小）时放弃此否定决策。

import numpy as np
import matplotlib.pyplot as plt


n = 5000            # Iterations (5000 years)
t = np.ones(n)
t[0] = 0            # Start time
dt = 1              # Time step

# Parameters
d0 = 200            # d0 (Start height)
s = 0.014           # Slope
lambd = 300         # lambda
sigma = 10000
xs = 40000          # x-coordinate of the "bump"
ac = 0.0005         # Accumulation
alfaf = 0.7         # Constant
alfam = 2           # Constant
eps = 1             # ...
delta = 1.127       # ...
c = 2.4             # Climate sensitivity

# Grids
x = np.linspace(0, 50000, 5000)     # My x grid
x[0] = 200                          # Start position of glacier's head
L = np.ones(n)                      # Glacier Length
B = np.zeros(n)
a = np.ones(n)*ac
Hf = np.ones(n)                     # Defining vectors to loop through
F = np.ones(n)
dLdt = np.zeros(n)

for i in range(1, n):
    if L[i - 1] < 0.0:
        # Since you are not allocating memory dynamically, you should discard the 
        # insignificant values of t and L (otherwise the chart will close).
        t = t[:i]
        L = L[:i]
        break
    t[i] = i
    a[i] = a[i - 1] + dt * ac                                                   # Accumulation (constant)
    x[i] = d0 - s * x[i] + lambd * np.exp(-((x[i] - xs) / sigma)**2)            # Bed geometry
    B[i] = a[i - 1] * L[i - 1]                                                  # Surface mass balance
    Hf[i] = np.max([alfaf * np.sqrt(L[i - 1]), -eps * delta * x[i - 1]])        # Ice thickness at glacier front
    F[i] = np.min([0, c * x[i - 1] * Hf[i - 1]])                                # Flux of ice at the glacier front
    dLdt[i] = (2 * (B[i - 1] + F[i - 1])) / (3 * alfam * np.sqrt(L[i - 1]))     # Differential equation
    L[i] = L[i - 1] + dt * dLdt[i - 1]                                          # Updating L with Eulers Method

plt.figure('Graph')
plt.plot(t, L)
plt.grid()
plt.show()