一次元波動方程式を、アニメーションにする方法

Question

pythonで、波動方程式の二次元関数を書いています。
z＝2(cos(x)) * (sin(cx))
c=1
とするとき、この二変数関数を二次元のプロットにしたいです・

コード自体にはエラーがないようなのですが、アニメーションにすると波が横に流れるものではなく、値の変わらない一本線になっていました。expの式で最初のテストを行った際は、右から左に波が移っていくアニメーションになっていたのですが、上記のcos,sinの形にするとうまくできませんでした。

#表示されるエラー
uvals.append(-2np.sin(t)np.sin(j))
uvals.append(U[t-1][j] + cdt/(2dx)*(U[t-1][j+1]-U[t-1][n-1]))

上の式を使うと、上下に動く式、下の式を使うと、直線が表示されるようになってしまいました。

initial の式を、
np.exp(-0.5*np.power(((x-0.5)/0.08), 2))
というものでプロットしたときは問題なく右から左にグラフが動く式になりました。動いたのですが、何故式を変えただけで

python

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# wave speed
c = 1

# spatial domain
xmin = 0
xmax = 1

n = 50 # num of grid points

# x grid of n points
X, dx = np.linspace(xmin,xmax,n,retstep=True)

# for CFL of 0.1
dt = 1/1024

# initial conditions
def initial_u(x):
    return 2*np.cos(x)*np.sin(0)


U = []


def u(x, t):
    if t == 0: # initial condition
        return initial_u(x)
    uvals = [] # u values for this time step
    for j in range(len(x)):
        if j == 0: # left boundary
            #上下に動く波になってしまいます
　　　　　　　#uvals.append(-2*np.sin(t)*np.sin(j))
　　　　　　　#ただの直線になってしまいます。
            #uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][n-1]))
        elif j == n-1: # right boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][0]-U[t-1][j-1]))
        else:　
　　　　　　　#上下に動く波になってしまいます
　　　　　　　#uvals.append(-2*np.sin(t)*np.sin(j))
　　　　　　　#ただの直線になってしまいます。
            #uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][j-1]))
    return uvals

# solve for 500 time steps
for t in range(500):
    U.append(u(X, t))

# plot solution
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)

# animate the time data
k = 0
def animate(i):
    global k
    x = U[k]
    k += 1
    ax1.clear()
    plt.plot(X,x)
    plt.grid(True)
    plt.ylim([-2,2])
    plt.xlim([0,1])

anim = animation.FuncAnimation(fig,animate,frames=360,interval=20)
plt.show()

#試してみたこと
下のものが、np.exp(-0.5*np.power(((x-0.5)/0.08), 2))
にて右から左に動くグラフです。
三角関数をプロットするためにいくつか変数は変えたのですが、何故下のプログラムは動くのに右のプログラムは動かないのかがわかりません。教えていただけないでしょうか。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# wave speed
c = 1

# spatial domain
xmin = 0
xmax = 1

n = 50 # num of grid points

# x grid of n points
X, dx = np.linspace(xmin,xmax,n,retstep=True)

# for CFL of 0.1
dt = 1/1024

# initial conditions
def initial_u(x):
    return np.exp(-0.5*np.power(((x-0.5)/0.08), 2))
    #return np.sin(x)*(1/10)

# each value of the U array contains the solution for all x values at each timestep
U = []

# explicit euler solution
def u(x, t):
    if t == 0: # initial condition
        return initial_u(x)
    uvals = [] # u values for this time step
    for j in range(len(x)):
        if j == 0: # left boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][n-1]))
        elif j == n-1: # right boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][0]-U[t-1][j-1]))
        else:
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][j-1]))
    return uvals

# solve for 500 time steps
for t in range(500):
    U.append(u(X, t))

# plot solution
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)

# animate the time data
k = 0
def animate(i):
    global k
    x = U[k]
    k += 1
    ax1.clear()
    plt.plot(X,x)
    plt.grid(True)
    plt.ylim([-2,2])
    plt.xlim([0,1])

anim = animation.FuncAnimation(fig,animate,frames=360,interval=20)
plt.show()

Comments

[hayataka2049]

コピペして実行してみましたが、
plt.plot(X,x)
のところで
ValueError: x and y must have same first dimension, but have shapes (50,) and (99,)
というエラーが発生しました。
とりあえず今一度コードをご確認ください。

[magichan]

hayataka2049さんの指摘の通り、x軸のデータ(X)サイズが50に対してy軸のデータ(x)のサイズが99となっているのでプロット時にエラーが出ます。
関数 u(x, t)　内にて j == n-1 以外の時に uvals.append() を２回ずつ行っておりますので、n=50 に対して uvals には 99 個のデータが格納されることが原因かと思われます。

[alpha_high]

>hayataka2049さま

コメントありがとうございます。
また、エラーについて教えて下さりありがとうございます。
修正し、無事グラフが出るようになりました。
私はこのグラフにて、グラフが右から左へ流れて行くものを想像していたのですが、現在は上下にのびたり縮んだりするグラフになってしまいました。
アルゴリズム 通りに書いたつもりだったのですが、私のコードでは、関数が動くという物にはなっていないのでしょうか。厚かましいお願いですが、見ていただけると本当にありがたいです。この度は、すぐにご回答頂き、本当にありがとうございました。

>magichanさま

コメントありがとうございます。また、私の問題の箇所について的確に教えて頂き、本当に感謝しています。本当にありがとうございます。無事、グラフは表示できるようになりました。しかし、私が当初予定していたもの（関数が右から左に動くだけのもの）が、上下に激しく動くグラフになってしまいました。波動方程式のアルゴリズムを使ったつもりだったのですが、もし何か私のコードで間違って居る点があれば、教えて頂きたいです。何度も申し訳ありませんが、教えて頂けると幸いです。
この度は、私のコードを見て頂き、本当にありがとうございました。

Answer 1

特に理屈はわかっていないんですが、2*np.cos(x)*np.sin(0)の初期値はぜんぶ0になるはずなので動きようがないのでは？

Comments

[alpha_high]

回答ありがとうございます。おっしゃる通りですね。初期値が０だと私のアルゴリズムだと動くはずがないです。試しに初期値が０にならない式で動かしてみたところ、きちんとグラフが表示されました。ありがとうございます。このアルゴリズムで、もし、初期値が０になるグラフを動かしたい場合は、初期値が０にならないような計算式に変換する、という方法意外にもありますでしょうか。他のアルゴリズムではなかなかきれいにシミュレーションできなかったので、このコードを使っていろんな計算式をアニメーションしたいのですが、厳しいでしょうか。もし宜しければ、どのような考え方をすれば、初期値が０でも表示できる、と教えて頂けると本当にありがたいです。

[hayataka2049]

初期値がすべて0というのはそもそも妥当な想定なのでしょうか。波動方程式に関してはそんなにしっかり勉強しなかったのと記憶があやふやなので、それ以上のことは言えませんが……