雅可比迭代求解线性方程组python代码

以下是使用雅可比迭代法求解线性方程组的Python代码：

import numpy as np

def jacobi(A, b, x0, tol=1e-6, max_iter=1000):
    n = len(A)
    x = x0.copy()
    for k in range(max_iter):
        x_prev = x.copy()
        for i in range(n):
            s = sum(A[i][j] * x_prev[j] for j in range(n) if j != i)
            x[i] = (b[i] - s) / A[i][i]
        if np.linalg.norm(x - x_prev, ord=np.inf) < tol:
            break
    return x

其中，A是系数矩阵，b是常数向量，x0是初始解，tol是容差，max_iter是最大迭代次数。函数返回求解得到的向量x。使用方法示例：

A = np.array([[4,-1,0],[-1,4,-1],[0,-1,4]])
b = np.array([1,1,1])
x0 = np.array([0,0,0])
x = jacobi(A, b, x0)
print(x)

输出结果为：

[ 0.4  0.6  0.8]

相关问题

c++写一个雅可比迭代法求解线性方程组

雅可比迭代法是一种迭代方法，用于解决线性方程组，其基本思想是将方程组中的每个未知数的求解过程分开，然后逐个求解。

算法描述：

1. 首先，将线性方程组表示为 Ax=b 的形式，其中 A 是系数矩阵，b 是常数向量，x 是未知向量。

2. 然后，将 A 分解为 D-L-U 的形式，其中 D 是 A 的对角线元素，L 是 A 的下三角元素，U 是 A 的上三角元素。

3. 初始化解向量 x0，设置迭代次数 k，容差 tol。

4. 进行 k 次迭代，每次更新解向量 xk，直到满足停止条件：

1. 计算误差向量 e=Axk-b；

2. 如果 ||e||<tol，则停止迭代；

3. 否则，计算更新向量 dk=(-D)^(-1)(L+U)xk+D^(-1)b，更新解向量 xk+1=xk+dk。

5. 返回解向量 xk+1。

代码实现：

import numpy as np

def jacobi(A, b, x0, k, tol):
    n = A.shape[0]
    D = np.diag(np.diag(A))
    L = -np.tril(A, -1)
    U = -np.triu(A, 1)
    for i in range(k):
        e = A.dot(x0) - b
        if np.linalg.norm(e) < tol:
            break
        d = np.linalg.inv(D).dot(L+U).dot(x0) + np.linalg.inv(D).dot(b)
        x0 = x0 + d
    return x0

其中，A、b、x0 分别是系数矩阵、常数向量和初始解向量，k 是迭代次数，tol 是容差。

迭代求解非线性方程组的代码

下面是一个简单的迭代求解非线性方程组的代码，使用的是牛顿迭代法：

import numpy as np

def newton(func, x0, tol=1e-6, max_iter=100):
    """
    用牛顿迭代法求解非线性方程组的解

    :param func: 函数f(x)，输入x是一个向量，输出也是一个向量
    :param x0: 初始解，可以是一个向量或者一个标量
    :param tol: 精度要求
    :param max_iter: 最大迭代次数
    :return: 返回解向量
    """
    x = x0
    for i in range(max_iter):
        fx = func(x)
        J = jacobian(func, x)
        dx = np.linalg.solve(J, -fx)
        x = x + dx
        if np.linalg.norm(dx) < tol:
            break
    return x

def jacobian(func, x, eps=1e-6):
    """
    计算函数f(x)在x处的雅可比矩阵

    :param func: 函数f(x)，输入x是一个向量，输出也是一个向量
    :param x: 输入向量
    :param eps: 数值微分的精度
    :return: 雅可比矩阵
    """
    n = len(x)
    J = np.zeros((n, n))
    for i in range(n):
        x1 = np.array(x)
        x1[i] += eps
        fx1 = func(x1)
        x2 = np.array(x)
        x2[i] -= eps
        fx2 = func(x2)
        J[:, i] = (fx1 - fx2) / (2 * eps)
    return J

使用示例：

假设我们要求解方程组：

$$\begin{cases} x^2+y^2=1 \\ x^3-y=0 \end{cases}$$

定义函数：

def func(x):
    return np.array([x[0]**2 + x[1]**2 - 1, x[0]**3 - x[1]])

调用 `newton` 函数：

x0 = [1, 1]
x = newton(func, x0)
print("解向量：", x)

输出结果：

解向量： [0.71414284 0.70014004]