daily pastebin goal
53%
SHARE
TWEET

Untitled

a guest Dec 11th, 2017 50 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import scipy.sparse as sp
  3. import matplotlib.pyplot as plt
  4.  
  5. def fun(x, t):
  6.     t_arr = np.zeros(len(x))
  7.     t_arr.fill(t)
  8.     return np.exp(-400 * ((x - t_arr)**2))
  9.  
  10. def adv(h, cfl):
  11.     num_x_vals = (3/h)
  12.     x = np.linspace(-1, 2, num_x_vals + 1)
  13.     delta_t = cfl * h
  14.     num_t_vals = (1/delta_t)
  15.     t = np.linspace(0, 1, num_t_vals + 1)
  16.     diagonals_center = np.zeros(len(x))
  17.     diagonals_bottom = np.zeros(len(x) - 1)
  18.     diagonals_top = np.zeros(len(x) - 1)
  19.     diagonals_bottom.fill(-1)
  20.     diagonals_top.fill(1)
  21.     tridiagonals = [diagonals_center, diagonals_bottom, diagonals_top]
  22.     C = sp.diags(tridiagonals, [0, -1, 1]).toarray()
  23.     C[len(C) - 1][0] = 1
  24.     C[0][len(C) - 1] = -1
  25.     C *= (-1/(2 * h))
  26.     f_0 = fun(x, 0)
  27.     f_1 = fun(x, 1)
  28.     f_neg1 = C @ fun(x, (-1 * delta_t))
  29.     f_neg2 = C @ fun(x, (-2 * delta_t))
  30.     u = np.zeros((len(x), len(t)))
  31.     u[:, 0] = f_0
  32.     u[:, 1] = u[:, 0] + (delta_t * (((23/12) * f_0) + ((-16/12) * f_neg1) + ((5/12) * f_neg2)))
  33.     for i in range(2, len(t)):
  34.         temp = f_neg1
  35.         f_neg2 = temp
  36.         f_neg1 = f_0
  37.         f_0 = C @ u[:, i - 1]
  38.         u[:, i] = u[:, i - 1] + (delta_t * (((23/12) * f_0) + ((-16/12) * f_neg1) + ((5/12) * f_neg2)))
  39.     err = np.amax(np.absolute(u[:, len(t) - 1] - f_1))
  40.     return f_1, u[:, len(t) - 1], x, err
  41.  
  42. u00, u01, x0, err0 = adv(0.02, 0.5)
  43. u10, u11, x1, err1 = adv(0.01, 0.5)
  44. u20, u21, x2, err2 = adv(0.005, 0.5)
  45. u30, u31, x3, err3 = adv(0.002, 0.5)
RAW Paste Data
Top