python-3.x - BicGStab 产生意外故障标志
问题描述
我需要解决一系列稀疏线性系统Ax=b
。第一个系统的解决方案x
是第二个系统的输入,第二个系统是第三个系统的输入,依此类推。由于数字误差复合和其他原因,我必须使用scipy.sparse.linalg.bicgstab
我的线性求解器。但是,对于一个甚至不是病态并且肯定有逆的系统,求解器会给我一个标志:“非法输入或故障”。
import numpy as np
from scipy.sparse.linalg import bicgstab, inv
from scipy import sparse
A = np.array(
[[ -1., 0., 0., 0., 0., 0., 0., 0.],
[ 0., -1., 0., 0., 0., 0., 0., 0.],
[ 0., 0., -10., 0., 0., 0., 0., 0.],
[ 0., 0., 0., -10., 0., 0., 0., 0.],
[ 0., 0., 3., 0., -3., 0., 0., 0.],
[ 0., 0., 0., 3., 0., -3., 0., 0.],
[ 0., 0., 0., 7., 3., 0., -10., 0.],
[ 0., 0., 7., 0., 0., 3., 0., -10.]]
)
A = -sparse.csc_matrix(A)
b = np.array([ 1., 0., 10., 0., 0., 0., 0., 0.])
x, flag = bicgstab(A=A, b=b, maxiter=40, tol=1e-6)
x, flag
>>> (array([1. , 0. , 1. , 0. , 1.00118012,
0. , 0.3004875 , 0.70009946]), -10)
只是为了证明这一点
inv(A).dot(b)
>>> array([1. , 0. , 1. , 0. , 1. , 0. , 0.3, 0.7])
上面的输出正是我所期望的。有谁知道为什么 bicgstab 没有给我想要的输出?我找不到有关 bicgstab 的非法输入或故障的文档,因此我是我的问题。
解决方案
-10
错误代码并不一定意味着您输入错误;在您的情况下,故障很可能是在迭代求解期间发生的。
通过稍微改变你的 RHS:
b = np.array([ 1., 0., 0., 0., 10., 0., 0., 0.])
scipy.bicgstab
即使没有前置条件,the也没有收敛的麻烦:
x, flag = bicgstab(A=A, b=b, maxiter=40, tol=1e-6)
print (x, flag)
(array([1. , 0. , 0. , 0. , 3.33333333, 0. , 1. , 0. ]), 0)
矩阵具有逆和合适的条件数的事实
print(np.linalg.cond(A))
14.616397823169317
不保证很容易获得特定 RHS的解决方案,尤其是使用迭代求解器或特定迭代求解器。在我看来(没有对矩阵谱及其内核空间进行详细分析),您的 RHS 恰好位于这样一个“糟糕的区域”。
如果您只是对解决方案感兴趣,我建议您切换到 GMRES:
x, flag = gmres(A=A, b=b, maxiter=40, tol=1e-6)
(数组([1. , 0. , 0.1 , 0. , 0.1 , 0. , 0.03, 0.07]), 0)
如果您有兴趣调查为什么 BiCGStab 失败,而 GMRES 在该系统的解决方案中取得成功,我会邀请您缩小范围的问题到Computational Science SE。
推荐阅读
- dynamics-crm - 自定义总金额
- c# - 如何将 TimeSpan 值的总和显示为 hh:mm:ss 格式的文本框
- java - 特定登录的 GET 方法
- keyboard - 将 win+n 映射到 X 和/或 i3 上的 Prior(PageUp)
- asp.net-web-api - 在每个 API 请求上重新生成 JWT 不记名令牌
- microsoft-graph-api - 如何使用 calendarView 查询事件冲突?
- asp.net-mvc - 部署到远程服务器 IIS 后如何访问 MVC 应用程序
- javascript - 跨域请求被阻止 - 服务器是否必须允许跨域请求?
- java - 将多个对象继承到一个类中
- c# - 您如何满足 WebApi 查询字符串中的 AND 和 OR 搜索?