Is it possible to use scikit-fem to simulate phi(x,t)_{tt}-psi(xxx)=0

[MisaelMil]

Is it possible to use scikit-fem to simulate this equation? If true, any help with writing the bilinear form is appreciated
$\begin{equation}\label{s1}
\rho_{1}\varphi(x,t){tt} - \kappa \psi(x,t){xxx}=0,
\end{equation} $
$(x,t) \in {(0,L) x (0,T)}.$
boundary and initial conditions can be defined as convenient

[kinnala]

Yes, I think it is possible. It looks like a wave equation with one extra partial derivative wrt x.

Perhaps you can just take this and replace L by something else: https://github.com/kinnala/scikit-fem/blob/8.1.0/docs/examples/ex44.py

Something like

def bilinf(u, v, w):
    return u.hess[0, 0] * v.grad[0]

and then use ElementLineHermite().

[MisaelMil]

Hi Thank you for your help. I tried the guidelines, however the graph did not represent the equation. Could you please see if I messed up any part of the code to insert the bilinear form? If the bilinear form is the same. if it's right then it must be numerical instability error. but didn't want to research the instability and be code error.

#beam

import numpy as np
from scipy.sparse import identity
from scipy.sparse.linalg import splu
from skfem import *
from skfem.models import laplace, mass
import matplotlib.pyplot as plt
import matplotlib.animation as anim
from pathlib import Path

dt = 0.01
theta = 1
nu=5
c = 1

m = MeshLine().refined(6)
basis = Basis(m, ElementLineHermite())

N = basis.N
L = laplace.assemble(basis)

#@BilinearForm
#def bilinf(u, v, w):
#from skfem.helpers import d,dd, dot, ddot
#return dot(d(u), d(v))
@BilinearForm
def bilinf(u, v, w):
return u.hess[0, 0] * v.grad[0]

M = asm(bilinf, basis)
I = identity(N)

A0 = bmat([[I, None],
[None, M]], "csc")

B0 = bmat([[None, I],
[-c ** 2 * L, None]], "csc")

A = A0 + theta * B0 * dt
B = A0 - (1. - theta) * B0 * dt
backsolve = splu(A).solve

def evolve(t, u):
while t < 1:
t = t + dt
u = backsolve(B @ u)
yield t, u

def bump(x):
#return np.exp(-100. * (x[0] - 0.5) ** 2)
return (np.sin(nunp.pix[0]))

U0 = basis.project(bump)
U = np.concatenate((U0, np.zeros(N)))

if name == "main":
fig, ax = plt.subplots()
ix = np.argsort(m.p[0])
line, = ax.plot(m.p[0, ix], U0[ix])

def update(t_u):
    t, u = t_u
    u1, _ = np.split(u, [M.shape[0]])
    line.set_ydata(u1[ix])

animation = anim.FuncAnimation(fig, update, evolve(0., U), interval=50, cache_frame_data=False)

output_file = Path(__file__).with_suffix(".gif")
#animation.save(output_file, writer='pillow')

plt.show()