Интегро-дифференциальное уравнение 1-го порядка для прогнозирования 3D-печати ⇐ Python

[Anonymous]

Прошу прощения, если это довольно простой вопрос. Я инженер и только начал изучать Python, поэтому, вероятно, упускаю что-то очень очевидное.
У меня есть интегро-дифференциальное уравнение, которое я пытаюсь решить с помощьюsolve_ivp для ОДУ и Quad для интеграла в ОДУ. Кажется, я не могу заставить его работать, я решил ОДУ с мю в качестве константы, но теперь мю является переменной, и вместо него нужно решать интегральное значение. Это уравнение, в котором mu является константой.


Константы:

Код: Выделить всё

kI_0 = 0.0015
μ = 3.74
z = np.linspace(0, 1.15, 100)

Код для решения этой проблемы:

Код: Выделить всё

import numpy as np
import math as ma
import matplotlib.pyplot as plt
import scipy as sp
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

#requires np.exp() overwise it doesnt work.
def dϕdt(t, ϕ, kI_0, μ, z):
return kI_0*(1-ϕ)*np.exp(-μ*z)

kI_0 = 0.0015
μ = 3.74

#Set the initial conditions as a vector of 100 zeros. Vector is required because its creating an array and isnt a scalar.
ϕ = np.zeros(100)

#linspace from 0 to 1.15 with 100 increment inbetween.
z = np.linspace(0, 1.15, 100)

#used t_span and dont define the increment size in t_span because then the function will select an appropriate size.
#use args to set the addional arguments with the ODE solver function.
sol_m1 = solve_ivp(dϕdt, y0=ϕ, t_span=(0, 3600), args=(kI_0, μ, z))

#Select the full row of the last column of y to get the finally results for
#y because I need the full conversion fraction value which is 1 when the full time has elapsed.
#this can be checked by changing the row (-1) to (1), this shows that only 1*10^-7 is the conversion fraction.
plt.plot(z, sol_m1.y[:,-1])
plt.ylabel("Conversion Fraction")
plt.xlabel("Moving Solid Liquid Interface (mm)")
plt.title("Liquid Conversion Fraction vs Depth \n (Numerically Solved)")
print(sol_m1)
#sol_m2

Этот код работает нормально и решает эту более простую ОДУ численно.
Однако, теперь, когда мю является переменной, уравнение выглядит следующим образом:
Изображение интеграло-дифференциала для решения
му можно записать в виде уравнения:

Код: Выделить всё

μ = μ_0*(1 - ϕ) + μ_infny*ϕ

Где φ получается из решенного ОДУ на каждом временном шаге.
Константы:

Код: Выделить всё

μ_0 = 1
μ_infny = 5
kI_0 = 0.0015

Я попытался написать следующий код:

Код: Выделить всё

import numpy as np
import math as ma
import matplotlib.pyplot as plt
import sympy as sy
import scipy as sp
from scipy.integrate import odeint
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from scipy.integrate import cumulative_trapezoid

μ_0 = 1
μ_infny = 5
kI_0 = 0.0015

def μ(ϕ, μ_0, μ_infny):
return μ_0*(1 - ϕ) + μ_infny*ϕ

#Set the initial conditions as a vector of 100 zeros. Vector is required because its creating an array and isnt a scalar.
ϕ = np.zeros(100)

#linspace from 0 to 1.15 with 100 increment inbetween.
z = np.linspace(0, 1.15, 100)

#create function for intergral

f = lambda z : -μ(ϕ, μ_0, μ_infny)
Int = quad(f, 0, 1.15)

#requires np.exp() overwise it doesnt work.
def dϕdt(t, ϕ, kI_0, Int):
return kI_0*(1-ϕ)*np.exp(Int)

#used t_span and dont define the increment size in t_span because then the function will select an appropriate size.
#use args to set the additional arguments with the ODE solver function.
sol_m1 = solve_ivp(dϕdt, y0=ϕ, t_span=(0, 3600), args=(kI_0, Int))

#Select the full row of the last column of y to get the finally results for
#y because I need the full conversion fraction value which is 1 when the full time has elapsed.
#this can be checked by changing the row (-1) to (1), this shows that only 1*10^-7 is the conversion fraction.
plt.plot(z, sol_m1.y[:,-1])
plt.ylabel("Conversion Fraction")
plt.xlabel("Moving Solid Liquid Interface (mm)")
plt.title("Liquid Conversion Fraction vs Depth \n (Numerically Solved)")
print(sol_m1)
#sol_m2

Есть ли у кого-нибудь идеи, как решить это численно?
Он должен создать график, подобный этому:
График того, что создает простой ODE
Из попытки кода это результат:
Ошибка типа: только массивы размера 1 могут быть преобразованы в скаляры Python