Toy examples for Nonstandard Finite Differences.

nonstandard_fd.png

nonstandard_fd.py

"""
Finite difference solution for the equation

  u' = -a*u
  u(0) = 1

We compare a standard finite difference with a non-standard one.

@author: Nicolás Guarín-Zapata
@date: May 2022

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

tmax = 1.0
a = 1.0
u0 = 1.0
npts_list = np.logspace(0.5, 5, 10, dtype=int)
h = tmax/npts_list
err_sfd = np.zeros_like(npts_list, dtype=float)
err_nsfd = np.zeros_like(npts_list, dtype=float)
for pt_cont, npts in enumerate(npts_list):
    dt = tmax/npts
    u_sfd = np.zeros((npts))
    u_nsfd = np.zeros((npts))
    u_sfd[0] = u_nsfd[0] = u0
    for cont in range(npts - 1):
        u_sfd[cont + 1] = u_sfd[cont] - a*dt*u_sfd[cont]
        u_nsfd[cont + 1] = u_nsfd[cont] - (1 - np.exp(-a*dt))*u_nsfd[cont]

    t = np.linspace(0, tmax, npts)
    u_e = np.exp(-a*t)
    err_sfd[pt_cont] = np.linalg.norm(u_e - u_sfd)
    err_nsfd[pt_cont] = np.linalg.norm(u_e - u_nsfd)


repo = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master"
style = repo + "/styles/minimalist.mplstyle"
plt.style.use(style)
plt.loglog(h, err_sfd, label="SFD")
plt.loglog(h, err_nsfd, label="NSFD")
plt.xlabel("$h$")
plt.ylabel("$\Vert u_e - u_h \Vert$")
plt.legend()
plt.savefig("nonstandard_fd.png", bbox_inches="tight", dpi=300)
plt.show()

nonstandard_fd_nl.png

nonstandard_fd_nl.py

"""
Finite difference solution for the equation

  u' - u² (1 - u) = cos(t) - sin²(t) + sin³(t)
  u(0) = 0

with solution

  u(t) = sin(t)

We compare a standard finite difference with a non-standard one.

@author: Nicolás Guarín-Zapata
@date: May 2022

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

tmax = 1
u0 = 0.0
f = lambda t: np.cos(t) - np.sin(t)**2 + np.sin(t)**3
npts_list = np.logspace(0.5, 5, 10, dtype=int)
h = tmax/npts_list
err_sfd = np.zeros_like(npts_list, dtype=float)
err_nsfd = np.zeros_like(npts_list, dtype=float)
for pt_cont, npts in enumerate(npts_list):
    dt = tmax/npts
    t = np.linspace(0, tmax, npts)
    u_e = np.sin(t)
    u_sfd = np.zeros((npts))
    u_nsfd = np.zeros((npts))
    u_sfd[0] = u_nsfd[0] = u0
    for cont in range(npts - 1):
        u_sfd[cont + 1] = u_sfd[cont] + dt*u_sfd[cont]**2*(1 - u_sfd[cont])\
                        + dt*f(t[cont])
        numer = u_nsfd[cont]*(1 + 2*(1 - np.exp(-dt))*u_nsfd[cont])\
            + (1 - np.exp(-dt)*f(t[cont]))
        denom = 1 + (1 - np.exp(-dt))*u_nsfd[cont]*(1 - u_nsfd[cont])
        u_nsfd[cont + 1] = numer/denom

    err_sfd[pt_cont] = np.linalg.norm(u_e - u_sfd)
    err_nsfd[pt_cont] = np.linalg.norm(u_e - u_nsfd)


repo = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master"
style = repo + "/styles/minimalist.mplstyle"
plt.style.use(style)
plt.loglog(h, err_sfd, label="SFD")
plt.loglog(h, err_nsfd, label="NSFD")
plt.xlabel("$h$")
plt.ylabel("$\Vert u_e - u_h \Vert$")
plt.legend()
plt.savefig("nonstandard_fd_nl.png", bbox_inches="tight", dpi=300)
plt.show()