2018-08-18-wave-eqns.ipynb

Kernel: SageMath 9.4

From Japanese translation of E. Kreyszig, Advanced Engineering Mathematics, 8th ed. （E. クライツィグ著，阿部寛治訳，「フーリエ解析と偏微分方程式」，培風館，p.168-170.） This is a computer simlatoin of a vibration of an unit sqare membrane.

A wave equation of 2 dimension $$\frac{\partial^2 u}{\partial t^2} = c^2 \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)$$ subject to boundary conditions $$u(x,y,t)=0 \text{ if $(x,y)ParseError: KaTeX parse error: Expected 'EOF', got '}' at position 20: …on the boundary}̲,\quad u(x,y,0)…$ the second and the third conditions are initial conditions.

Now we use a method of 'separation of variables'. Let $F(x,y)$ and $G(t)$ be a spatial part and a time dependent part of the solutoin $u(x,y,t)$: $u(x, y, t) = F(x,y)G(t)$. By a standard procedure, we have an ODE $\ddot{G}+\lambda G=0$ (time dependent part), and a PDE $F_{{xx}} + F_{yy} + \nu^{2}F=0$ (spatial part). The dot notation means that a differential with respect to $t$, and the subscript means a partial derivative with respect to $x$ or $y$.

The time dependent part is governed by the eigenvalue $$\lambda_{m,n} = c\pi \sqrt{m^{2}+n^{2}}.$$ For a $\lambda_{m,n}$, we have an eigenfunctoin $$G_{m,n}(t) = B_{m,n}\cos(\lambda_{m,n}t) + B_{m,n}^{*}\sin(\lambda_{m,n}t).$$ So the solution of the original equation which corresponds to $\lambda_{m,n}$ is $$u_{m,n}(x,y,t) = (B_{m,n}\cos(\lambda_{m,n}t) + B_{m,n}^{*}\sin(\lambda_{m,n}t)) \sin(m\pi x)\sin(n\pi y).$$

The general solution of the equation is given as a superposition of $u_{m,n}$. $$u(x,y,t) = \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} u_{m,n}(x,y,t) = \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} (B_{m,n}\cos(\lambda_{m,n}t) + B_{m,n}^{*}\sin(\lambda_{m,n}t)) \sin(m\pi x)\sin(n\pi y).$$ The coefficients $B_{m,n}$ and $B_{m,n}^{*}$ is given by the integrals $$B_{m,n} = 4\int_{0}^{1}\int_{0}^{1} f(x,y)\sin(m\pi x)\sin(n\pi y)\,dx\,dy$$ and $$B_{m,n}^{*} = \frac{4}{\lambda_{m,n}}\int_{0}^{1}\int_{0}^{1} g(x,y)\sin(m\pi x)\sin(n\pi y)\,dx\,dy.$$

x,y = var('x,y')

def b(m,n,f):
    """The coefficient of mode (m, n)."""
    return(4*numerical_integral(numerical_integral(f(x, y)*sin(m*pi.n()*x), 0, 1, params=[1])[0]*sin(n*pi.n()*y),0, 1)[0])

b(5,8, lambda x, y: (4*x-x^2)*(2*y-y^2))

-6.377494332652553e-17

def u(t, m, n, f):
    """The mode (m, n) eigen function."""
    return(b(m, n, f)*cos(pi.n()*sqrt(m^2 + n^2).n()*t)*sin(m*pi.n()*x)*sin(n*pi.n()*y))

u(0.1, 5, 13, lambda x, y: (4*x-x^2)*(2*y-y^2))

-0.0124248487328100*sin(15.7079632679490*x)*sin(40.8407044966673*y)

A nodal curve of an eigenfunction of mode $(4,7)$.

iplt=implicit_plot(u(0.1, 4, 7, lambda x, y: (4*x-x^2)*(2*y-y^2)), (x, 0.0, 1.0), (y, 0.0, 1.0)); iplt.show()

A density plot of the same function.

dplt= density_plot(u(0.1, 4, 7, lambda x, y: (4*x-x^2)*(2*y-y^2)), (x, 0.0, 1.0), (y, 0.0, 1.0), cmap='jet'); dplt.show()

(iplt+dplt).show()

It is interesting that the same eigenvalue $\lambda_{m,n}$ is attained by some different pairs $(m,n)$ of integers. For example, $\lambda_{4,7} = \lambda_{1,8}$ holds since $4^2+7^2 = 65 = 1^2+8^2$. This happens when $m^2+n^2$ have two or more prime factors $p$ and $q$ which are congruent to $1$ modulo $4$. These pairs of integers provide different eigenfunctions $u_{m,n}(t)$ associated to the same eigenvalue $\lambda_{m,n}$. The followings are examples of such eigenfunctions.

gr = plot3d(u(0.1, 4, 7, lambda x, y: (4*x-x^2)*(2*y-y^2)), (x, 0.0, 1.0), (y, 0.0, 1.0), plot_points=80)

gr.show(viewer='canvas3d')

MIME type unknown not supported

gr2 = plot3d(u(0.1, 5, 13, lambda x, y: (4*x-x^2)*(2*y-y^2)), (x, 0.0, 1.0), (y, 0.0, 1.0), plot_points=80)

gr2.show(viewer='jmol')

MIME type unknown not supported

def w(t, f):
    """The superposition u(t, f) of u(t, m, n, f)."""
    return(sum(sum(u(t, m, n, f) for m in range(10)) for n in range(10)))

gr = plot3d(w(0.1,  lambda x, y: (4*x-x^2)*(2*y-y^2)), (x, 0, 1), (y, 0, 1), plot_points=80)

gr.show(viewer='jmol', aspect_ratio=[1,1,0.3])

MIME type unknown not supported

To make an animation, first we prepare frames of animation.

def frames(t, m, n):
    return(plot3d(u(t, m, n, lambda x, y: (4*x-x^2)*(2*y-y^2))*10^17, (x, 0, 1), (y, 0, 1)))
    

I'm not sure that the following method is the best way to restrict a region to be drawn. (To use implicit plot).

z = var('z')
def framesw(t):
    return(implicit_plot3d(w(t, lambda x, y: (4*x-x^2)*(2*y-y^2))-z, (x, 0.0, 1.0), (y, 0.0, 1.0), (z, -6.0, 6.0), plot_points=80))

fs= [frames(t, 2, 3) for t in sxrange(0.0, 2.0, 0.1)]

animate(fs)

fsw = [framesw(t) for t in sxrange(0.0, 4.0, 0.05)]

anmw = animate(fsw)

anmw.show() # ??? does not work on cocalc.

WARNING: Some output was deleted.

anmw.apng("99wave0819.apng")

anmw.save("99wave0819.gif")

---

There is a different way to make an animation. I generate many pictures, then I make an animation by using 'convert (in ImageMagick)' on the command line.

gr = sage.plot.plot3d.parametric_surface.ParametricSurface(); gr
c = 0
for t in srange(0.0, 2.0, 0.1):
    gr = implicit_plot3d(u(t, 1.0, 2.0, lambda x,y: (4*x-x^2)*(2*y-y^2))-z, (x, 0.0, 1.0), (y, 0.0, 1.0), (z, -1.0e-16, 1.0e-16), opacity=0.8, adoptive=True)
    gr.save_image("99wave%s.png"%t)
print("done!")

done!