CoCalc Public Files2018-08-18-wave-eqns.ipynbOpen with one click!
Author: Iwao KIMURA
Views : 113
Description: a simulation of vibrating square membrane.
Compute Environment: Ubuntu 18.04 (Deprecated)

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 2ut2=c2(2ux2+2uy2)\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 if (x,y) is on the boundary,u(x,y,0)=f(x,y),utt=0=g(x,y),u(x,y,t)=0 \text{ if $(x,y)$ is on the boundary},\quad u(x,y,0)=f(x,y), \quad \left.\frac{\partial u}{\partial t} \right|_{t=0} = g(x,y), the second and the third conditions are initial conditions.

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

The time dependent part is governed by the eigenvalue λm,n=cπm2+n2. \lambda_{m,n} = c\pi \sqrt{m^{2}+n^{2}}. For a λm,n\lambda_{m,n}, we have an eigenfunctoin Gm,n(t)=Bm,ncos(λm,nt)+Bm,nsin(λm,nt). 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 λm,n\lambda_{m,n} is um,n(x,y,t)=(Bm,ncos(λm,nt)+Bm,nsin(λm,nt))sin(mπx)sin(nπy). 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 um,nu_{m,n}. u(x,y,t)=m=1n=1um,n(x,y,t)=m=1n=1(Bm,ncos(λm,nt)+Bm,nsin(λm,nt))sin(mπx)sin(nπy). 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 Bm,nB_{m,n} and Bm,nB_{m,n}^{*} is given by the integrals Bm,n=40101f(x,y)sin(mπx)sin(nπy)dxdy B_{m,n} = 4\int_{0}^{1}\int_{0}^{1} f(x,y)\sin(m\pi x)\sin(n\pi y)\,dx\,dy and Bm,n=4λm,n0101g(x,y)sin(mπx)sin(nπy)dxdy. 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.

In [1]:
x,y = var('x,y')
In [2]:
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])
In [3]:
b(5,8, lambda x, y: (4*x-x^2)*(2*y-y^2))
-6.377494332652553e-17
In [4]:
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))
In [5]:
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)(4,7).

In [6]:
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.

In [7]:
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()
In [8]:
(iplt+dplt).show()

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

In [9]:
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)
In [10]:
gr.show(viewer='canvas3d')
In [11]:
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)
In [12]:
gr2.show(viewer='jmol')