Пока l oop для радиоактивного распада. ValueError: установка элемента массива с последовательностью

Question

def calculateEuler(nuclei_A, nuclei_B, tau_A, tau_B, time, dt, n):
    i = 0
    while i < n-1:
        nuclei_A[i+1] = nuclei_A[i] - (nuclei_A[i]/tau_A)*dt
        nuclei_B[i+1] = nuclei_B[i] + (nuclei_A[i]/tau_A)*dt - (nuclei_B[i]/tau_B)*dt
        time[i+1] = time[i] + dt
        i = i+1
    return nuclei_A, nuclei_B, time

Я пытаюсь написать функцию, которая аппроксимирует количество ядер (А и В), когда они подвергаются go радиоактивному распаду. Начальное число ядер, nuclei_A [0] и nuclei_B [0] уже было назначено в предыдущей функции.

Я получаю две ошибки для линии, определяющей nuclei_B.:

TypeError: only size-1 arrays can be converted to Python scalars.

ValueError: setting an array element with a sequence.

Мне трудно понять, что означают эти ошибки и что я делаю неправильно.

Answer 1

Следующее проходит без ошибок, и каждый из видов достигает нуля за отведенное время. Я не внес никаких изменений в вашу функцию. Если вы обнаружите, что инициализация с другими значениями нарушает этот код, вы можете переформатировать свой вопрос с помощью этого примера кода и этих значений.

import numpy as np
import matplotlib.pyplot as plt

def calculateEuler(nuclei_A, nuclei_B, tau_A, tau_B, time, dt, n):
    i = 0
    while i < n-1:
        nuclei_A[i+1] = nuclei_A[i] - (nuclei_A[i]/tau_A)*dt
        nuclei_B[i+1] = nuclei_B[i] + (nuclei_A[i]/tau_A)*dt - (nuclei_B[i]/tau_B)*dt
        time[i+1] = time[i] + dt
        i = i+1
    return nuclei_A, nuclei_B, time

n = 100
A = np.zeros((n))
B = np.zeros((n))
A[0] = 10000
B[0] = 10000
tA = 1.6
tB = 2
time = np.zeros((n))
dt = .5

rA, rB, rT = calculateEuler(A, B, tA, tB, time, dt, n)

plt.figure()
plt.plot(rT, rA)

plt.figure()
plt.plot(rT, rB)

Я немного повеселился с этим и смог существенно ускорить код только по простой причине. Вышеописанная функция выполняется в среднем за ~ 2,91 мс для имитации с шагом 1000 раз. Приведенная ниже функция выполняется в среднем за ~ 0,009 мс для того же массива времени.

import numpy as np
import matplotlib.pyplot as plt
from numba import njit

@njit(nogil=True, cache=True)
def calculateEuler(nuclei_A, nuclei_B, tau_A, tau_B, dt):

    for idx in range(nuclei_A.size):
        if idx == nuclei_A.size - 1 : break
        nuclei_A[idx+1] = nuclei_A[idx] - (nuclei_A[idx]/tau_A)*dt
        nuclei_B[idx+1] = nuclei_B[idx] + (nuclei_A[idx]/tau_A)*dt - (nuclei_B[idx]/tau_B)*dt
    return nuclei_A, nuclei_B

n = 1000
A = np.zeros((n))
B = np.zeros((n))
A[0] = 1e7
B[0] = 1e7
tA = 1.6
tB = 2
dt = 0.5
time = np.arange(0, n*dt, dt)

rA, rB = calculateEuler(A, B, tA, tB, dt)