Одномерное волновое уравнение - PullRequest
Новая
       62

Одномерное волновое уравнение

0 голосов
/ 24 апреля 2020

Я пытался понять, как решить одномерное волновое уравнение, используя метод конечных разностей и интегрируя с Рунге-Кутта 4-го порядка (не с помощью функции). Начальные условия: u0 (x) = 1, u0 (xcenter) = 10. Не могли бы вы помочь мне с этим, пожалуйста?

enter image description here

Я смотрел на некоторые сценарии, но не понимают все части. Например: 

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

dx=0.1 #space increment
dt=0.05 #time increment
tmin=0.0 #initial time
tmax=3.0 #simulate until
xmin=-5.0 #left bound
xmax=5.0 #right bound...assume packet never reaches boundary
c=1.0 #speed of sound
rsq=(c*dt/dx)**2 #appears in finite diff sol - pravá strana

nx = int((xmax-xmin)/dx) + 1 #number of points on x grid
nt = int((tmax-tmin)/dt) + 2 #number of points on t grid
u = np.zeros((nt,nx)) #solution to WE

#set initial pulse shape
def init_fn(x):
    val = 10
    if val<.001:
        return 1
    else:
        return val

for a in range(0,nx):
    u[0,a]=init_fn(xmin+a*dx)
    u[1,a]=u[0,a]

#simulate dynamics
for t in range(1,nt-1):
    for a in range(1,nx-1):
        u[t+1,a] = 2*(1-rsq)*u[t,a]-u[t-1,a]+rsq*(u[t,a-1]+u[t,a+1])

fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
for i in range(nt):
    p, = plt.plot(u[i,:], 'k')   # this is how you'd plot a single line...
    plts.append( [p] )           # ... but save the line artist for the animation
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat_delay=300)   # run the animation
ani.save('wave.mp4')    # optionally save it to a file

plt.show()

Я потерян из функции def init_fn (x): Это означает начальные условия?

Здесь я не знаю, как определить x и t. Я попытался сделать это так же, как и в предыдущем сценарии, но он не работает. 

# Given mesh points as arrays x and t (x[i], t[n])
dx = x[1] - x[0]
dt = t[1] - t[0]
C = c*dt/dx            # Courant number
Nt = len(t)-1
C2 = C**2              # Help variable in the scheme

# Set initial condition u(x,0) = I(x)
for i in range(0, Nx+1):
    u_1[i] = I(x[i])

# Apply special formula for first step, incorporating du/dt=0
for i in range(1, Nx):
    u[i] = u_1[i] + 0.5*C**2(u_1[i+1] - 2*u_1[i] + u_1[i-1])
u[0] = 0;  u[Nx] = 0   # Enforce boundary conditions

# Switch variables before next step
u_2[:], u_1[:] = u_1, u

for n in range(1, Nt):
    # Update all inner mesh points at time t[n+1]
    for i in range(1, Nx):
        u[i] = 2*u_1[i] - u_2[i] + C**2(u_1[i+1] - 2*u_1[i] + u_1[i-1])

    # Insert boundary conditions
    u[0] = 0;  u[Nx] = 0

    # Switch variables before next step
    u_2[:], u_1[:] = u_1, u

Здесь я попытался написать, какие части он мог бы иметь, но я не могу закончить sh это. 

import numpy as np 
from numpy.linalg import inv
from matplotlib import pyplot as plt

# -------------------------------------------------------------------
# initial conditions
# -------------------------------------------------------------------

F = np.array([0.0,0.0])
print(F)

def F(t):
    F = np.array([0.0,0.0])
    return F

def G(y,t): return A_inv.dot( F(t) - B.dot(y) )


# -------------------------------------------------------------------
# integrator
# -------------------------------------------------------------------

def RK4_step(y, t, dt):
    k1 = G(y,t)
    k2 = G(y+0.5*k1*dt, t+0.5*dt)
    k3 = G(y+0.5*k2*dt, t+0.5*dt)
    k4 = G(y+k3*dt, t+dt)

    return dt * (k1 + 2*k2 + 2*k3 + k4) /6

# variables
m = 2.0
k = 2.0
c = 0.0   # critical damping = 2 * SQRT(m*k) = 4.0

F0 = 0.0
delta_t = 0.1
omega = 1.0
time = np.arange(0.0, 40.0, delta_t)

# initial state
y = np.array([0,1])   # [velocity, displacement]

A = np.array([[m,0],[0,1]])
B = np.array([[c,k],[-1,0]])
A_inv = inv(A)

Y = []
force = []

# time-stepping solution
for t in time:
    y = y + RK4_step(y, t, delta_t) 

    Y.append(y[1])
    force.append(F(t)[0])

    KE = 0.5 * m * y[0]**2
    PE = 0.5 * k * y[1]**2

    if t % 1 <= 0.01:
        print 'Total Energy:', KE+PE

# plot the result
plt.plot(time,Y)
plt.plot(time,force)
plt.grid(True)
plt.legend(['Displacement', 'Force'], loc='lower right')
plt.show()



Lx = 10
dx = 0.1
# odhad kroku
nx = np.fix(Lx/dx)
x = linspace(0, Lx, nx)

# Time 
T = 35

# Tři proměnné - v čase n-1 a v čase n+1
wn = zeros(nx,1)
wnm1 = wn
wnp1 = wn

# dt value
CLF = 1
c = 1
dt = CFL*dx/c

# Initial conditions
wn[1/2*nx-1:1/2*nx+1] = 10
wp1[:] = wn(:);


# Boundary conditions
# Kroky dle času
t = 0      

while(t < T):
   wn[0] = 0
   wn[-1] = 0


# Solution
t = t+dt
wnm1 = wn
wn = wnp1

Большое спасибо

...