У меня есть следующий код. Начало довольно долгое, но служит только для генерации данных. Проблема возникает с несколькими строками в конце.
##### Import packages
import numpy as np
import scipy.linalg as la
##### Initial conditions
N = 5
lamda = 7
mu = 2
a = 0.5
r = - np.log(a).copy()
St_Sp = np.arange(- N, N + 1)
Card = St_Sp.shape[0]
##### Define infintesimal generator
def LL(x, y):
if x == N or x == - N: re = 0
elif x - y == - 1: re = lamda
elif x - y == 1: re = mu
elif x - y == 0: re = - (mu + lamda)
else: re = 0
return re
def L(x):
return - LL(x, x)
##### Define function Phi
def Phi(x): return max(x, 0)
Phi = np.vectorize(Phi)
##### Define vector b
b = Phi(St_Sp).copy()
##### Define function Psi
def Psi(x): return L(x) / (L(x) + r)
Psi = np.vectorize(Psi)
##### Generate a Boolean vector whose all elements are False
d = np.array([0] * Card).astype(bool)
##### Define matrix A
A = np.zeros((Card, Card))
for i in range(Card):
for j in range(Card):
if (i != j) & (L(St_Sp[i]) != 0):
A[i, j] = LL(St_Sp[i], St_Sp[j]) / L(St_Sp[i])
elif (i != j) & (L(St_Sp[i]) == 0):
A[i, j] = 0
elif (i == j) & (Psi(St_Sp[i]) != 0):
A[i, j] = - 1 / Psi(St_Sp[i])
else: A[i, j] = 1
##### Row names of A
rows = np.arange(0, Card)
##### Define matrix B
B = np.zeros((Card, Card))
for i in range(Card):
for j in range(Card):
if i != j:
B[i, j] = LL(St_Sp[i], St_Sp[j])
else: B[i, j] = LL(St_Sp[i], St_Sp[j]) - r
##### Generate I_0
I = [np.array([1] * Card).astype(bool), d.copy()]
Z = b.copy()
index0 = np.matmul(B, Z) <= 0
index1 = ~ index0
##### Generate I_1
I = [index0, index1]
Z = b.copy()
if np.sum(I[1]) > 0:
order = np.concatenate((rows[I[1]], rows[~ I[1]]))
A1 = A[np.ix_(rows[I[1]], order)]
A2 = la.lu(A1)[2]
p = np.atleast_2d(A1).shape[0]
B1 = A2[:, range(p)]
B2 = - np.matmul(A2[:, p:], Z[I[0]])
print('Before being assigned new values, Z is \n', Z)
print('\n The index I[1] of elements of Z to be change \n', I[1])
M = la.solve_triangular(B1, B2, lower = False)
print('\n The values to be assigned to Z[I[1]] is \n', M)
Z[I[1]] = M
print('\n After being assigned new values, Z is \n', Z)
с результатом
Before being assigned new values, Z is
[0 0 0 0 0 0 1 2 3 4 5]
The index I[1] of elements of Z to be change
[False False False False False True True True True True False]
The values to be assigned to Z[I[1]] is
[2.08686055 2.88974949 3.40529229 3.88978577 4.41338306]
After being assigned new values, Z is
[0 0 0 0 0 2 2 3 3 4 5]
Мне очень странно, что команда Z[I[1]] = M не назначает новые значения из M на позицию Z с индексом I[1]. Не могли бы вы рассказать, почему возникает эта проблема и как ее решить?