Python linalg.spsolve方法代码示例

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def test_set_up_tms(self, tms_sphere):
        m, cond, dAdt, E_analytical = tms_sphere
        S = fem.FEMSystem.tms(m, cond)
        b = S.assemble_tms_rhs(dAdt)
        x = spalg.spsolve(S.A, b)
        v = mesh_io.NodeData(x, 'FEM', mesh=m)
        E = -v.gradient().value * 1e3 - dAdt.node_data2elm_data().value
        m.elmdata = [mesh_io.ElementData(E_analytical, 'analytical'),
                     mesh_io.ElementData(E, 'E_FEM'),
                     mesh_io.ElementData(E_analytical + dAdt.node_data2elm_data().value, 'grad_analytical'),
                     dAdt]

        m.nodedata = [mesh_io.NodeData(x, 'FEM')]
        #mesh_io.write_msh(m, '~/Tests/fem.msh')
        assert rdm(E, E_analytical) < .2
        assert np.abs(mag(E, E_analytical)) < np.log(1.1) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def test_dirichlet_problem_cube(self, cube_lr):
        m = cube_lr.crop_mesh(5)
        cond = np.ones(len(m.elm.tetrahedra))
        top = m.nodes.node_number[m.nodes.node_coord[:, 2] > 49]
        bottom = m.nodes.node_number[m.nodes.node_coord[:, 2] < -49]
        bcs = [fem.DirichletBC(top, np.ones(len(top))),
               fem.DirichletBC(bottom, -np.ones(len(bottom)))]
        S = fem.FEMSystem(m, cond)
        A = S.A
        b = np.zeros(m.nodes.nr)
        dof_map = S.dof_map
        for bc in bcs:
            A, b, dof_map = bc.apply(A, b, dof_map)
        x = spalg.spsolve(A, b)
        for bc in bcs:
            x, dof_map = bc.apply_to_solution(x, dof_map)
        order = dof_map.inverse.argsort()
        x = x[order]
        sol = m.nodes.node_coord[:, 2]/50.
        assert np.allclose(sol, x.T) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve ( self, A, b ):

    if len(A.shape) == 2:
      C = self.getConstraintsMatrix()

      a = zeros(len(self))
      a[self.constrainedDofs] = self.constrainedFac * array(self.constrainedVals)

      A_constrained = C.transpose() * (A * C )
      b_constrained = C.transpose()* ( b - A * a )

      x_constrained = spsolve( A_constrained, b_constrained )

      x = C * x_constrained

      x[self.constrainedDofs] = self.constrainedFac * array(self.constrainedVals)
    
    elif len(A.shape) == 1:
      x = b / A

      x[self.constrainedDofs] = self.constrainedFac * array(self.constrainedVals)
   
    return x 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve(self):

        shape = self.bcs.shape
        if self._L is None:
            self._L = self.lhs.matrix(shape) # expensive operation, so cache it

        L = sparse.lil_matrix(self._L)
        f = self.rhs.reshape(-1, 1)

        nz = list(self.bcs.row_inds())

        L[nz, :] = self.bcs.lhs[nz, :]
        f[nz] = np.array(self.bcs.rhs[nz].toarray()).reshape(-1, 1)

        L = sparse.csr_matrix(L)
        return spsolve(L, f).reshape(shape) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = mesh(self.nvars)
        me.values = LA.spsolve(sp.eye(self.nvars)-factor*self.A,rhs.values)

        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """
        
        M  = self.Id - factor*self.A
        
        b = np.concatenate( (rhs.values[0,:], rhs.values[1,:]) )
        
        sol = LA.spsolve(M, b)

        me = mesh(self.nvars)
        me.values[0,:], me.values[1,:] = np.split(sol, 2)
        
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = mesh(self.nvars)
        me.values = LA.spsolve(sp.eye(self.nvars)-factor*self.A,rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system_jacobian(self, dfdu, rhs, factor, u0, t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            dfdu: the Jacobian of the RHS of the ODE
            rhs: right-hand side for the linear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = self.dtype_u(self.init)
        me.values = spsolve(sp.eye(self.params.nvars) - factor * dfdu, rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)

        if self.params.direct_solver:
            me.values = spsolve(self.Id - factor * self.A, rhs.values.flatten())
        else:
            me.values = gmres(self.Id - factor * self.A, rhs.values.flatten(), x0=u0.values.flatten(),
                              tol=self.params.lintol, maxiter=self.params.liniter)[0]
        me.values = me.values.reshape(self.params.nvars)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def timestep(self, u0, dt):

        # Solve for stages
        for i in range(0, self.nstages):

            # Construct RHS
            rhs = np.copy(u0)
            for j in range(0, i):
                rhs += dt * self.A_hat[i, j] * (self.f_slow(self.stages[j, :])) + dt * self.A[i, j] * \
                    (self.f_fast(self.stages[j, :]))

            # Solve for stage i
            if self.A[i, i] == 0:
                # Avoid call to spsolve with identity matrix
                self.stages[i, :] = np.copy(rhs)
            else:
                self.stages[i, :] = self.f_fast_solve(rhs, dt * self.A[i, i])

        # Update
        for i in range(0, self.nstages):
            u0 += dt * self.b_hat[i] * (self.f_slow(self.stages[i, :])) + dt * self.b[i] * \
                (self.f_fast(self.stages[i, :]))

        return u0 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the nonlinear system
            factor (float): abbrev. for the node-to-node stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver (not used here so far)
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        M = self.Id - factor * self.A

        b = np.concatenate((rhs.values[0, :], rhs.values[1, :]))

        sol = spsolve(M, b)

        me = self.dtype_u(self.init)
        me.values[0, :], me.values[1, :] = np.split(sol, 2)

        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        self.uext.values[0] = 0.0
        self.uext.values[-1] = 0.0
        self.uext.values[1:-1] = rhs.values[:]
        me.values = spsolve(sp.eye(self.params.nvars + 2, format='csc') - factor * self.A, self.uext.values)[1:-1]
        # me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A[1:-1, 1:-1], rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system_1(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(u0)
        me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A, rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A, rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import spsolve [as 别名]
def solve_direct_lagrange(k_lg, f):
    """Solves a linear system of equations (Ku = f) using the direct solver method and the
    Lagrangian multiplier method.

    :param k: (N+1) x (N+1) Lagrangian multiplier matrix of the linear system
    :type k: :class:`scipy.sparse.csc_matrix`
    :param f: N x 1 right hand side of the linear system
    :type f: :class:`numpy.ndarray`

    :return: The solution vector to the linear system of equations
    :rtype: :class:`numpy.ndarray`

    :raises RuntimeError: If the Lagrangian multiplier method exceeds a tolerance of 1e-5
    """

    u = spsolve(k_lg, np.append(f, 0))

    # compute error
    err = u[-1] / max(np.absolute(u))

    if err > 1e-5:
        err = "Lagrangian multiplier method error exceeds tolerance of 1e-5."
        raise RuntimeError(err)

    return u[:-1]