Heat equation: I can't find this method "penalize" is not defined

[MisaelMil]

r"""Heat equation in one-dimension.

The solution of the heat equation is reduced from two dimensions in ex19 to
one here, to demonstrate the different post-processing required.

The two-dimensional modes from ex19 reduce in the limit
:math:w_1 \rightarrow\infty of the strip :math:|x| < w to

.. math::
\exp \left[
-\kappa
\left{
\left(\frac{(2n + 1)\pi}{2w}\right)^2
\right}
t
\right]
\cos\frac{(2n + 1)\pi x}{2w}
for :math:n = 0, 1, 2, \ldots.

Here we simulate the decay of the fundamental, :math:n = 0, again
discretizing time using the generalized ('theta method') trapezoidal
rule.

"""
from typing import Iterator, Tuple

import numpy as np
from scipy.sparse.linalg import splu

from skfem import *
from skfem.models.poisson import laplace, mass

import matplotlib.pyplot as plt

#imputs____________________________________
halfwidth = 2.0
ncells = 2 ** 3
diffusivity = 5.0
dt = 0.01
theta = 0.5 # Crank–Nicolson
#__________________________________________

mesh = MeshLine(np.linspace(-1, 1, 2 * ncells) * halfwidth)

element = ElementLineP2() # or ElementLineP1
basis = Basis(mesh, element)

L = diffusivity * asm(laplace, basis)
M = asm(mass, basis)

print("dt =", dt)

L0, M0 = penalize(L, M, D=basis.get_dofs())
A = M0 + theta * L0 * dt
B = M0 - (1 - theta) * L0 * dt

backsolve = splu(A.T).solve # .T as splu prefers CSC

u_init = np.cos(np.pi * basis.doflocs[0] / 2 / halfwidth)

def exact(x: float, t: float) -> float:
return np.exp(-diffusivity * (np.pi / 2 / halfwidth) ** 2 * t) * np.cos(np.pi * x / 2 / halfwidth)

def evolve(t: float, u: np.ndarray) -> Iterator[Tuple[float, np.ndarray]]:
while np.linalg.norm(u, np.inf) > 2 ** -3:
t, u = t + dt, backsolve(B @ u)
yield t, u

#post processing____________________________
if name == "main":
from mpl_toolkits.mplot3d import Axes3D

# Initialize an array to store numerical solution values
num_time_steps = 50
surface_matrix = np.zeros((num_time_steps, len(basis.doflocs[0])))

# Calculate numerical solution for different time points and x positions
for i, (t, u) in enumerate(evolve(0, u_init)):
    if i < num_time_steps:
        surface_matrix[i, :] = u

#Create the X and Y Axes
X, Y = np.meshgrid(basis.doflocs[0], np.linspace(0, 1, num_time_steps))

# Create the 3D Surface Graphic
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, surface_matrix, cmap='viridis')

ax.set_xlabel('x')
ax.set_ylabel('Time')
ax.set_zlabel('numerical solution')
plt.show()

[kinnala]

Here it is:

scikit-fem/skfem/utils.py

Line 403 in 760786c

def penalize(A: spmatrix,

[MisaelMil]

Thank you ! :) Enviado do Yahoo Mail no Android Em sáb., 26 26e ago. 26e 2023 às 17:33, ***@***.***> escreveu: Closed #1052 as resolved. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>