Finite Differences

The Nature of Mathematical Modeling (draft)

ODEs

  • Taylor series as numerical method
  • order of expansion
  • Euler's method
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .6
width = 10*scale
height = 5*scale
marker = 2*height
linew = marker/2.5
font = 5.5*height
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def line(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=0.8*linew)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/7,headwidth=marker,headlength=marker))
arrow(1,1,1,4)
arrow(1,1,8,1)
line(2.5,.8,2.5,1.2)
line(6.5,.8,6.5,1.2)
point(2.5,2.5)
point(6.5,2.01)
point(6.5,3.5)
text(2.5,0.5,'$x$')
text(6.5,0.5,'$x+h$')
text(1,4.35,'$y$')
arrow(2.45,2.47,6.35,3.45)
def curve(x0,x1,a,b,c):
    x0 = scale*x0
    x1 = scale*x1
    x = np.linspace(x0,x1,100)
    y = a+b*x+c*x*x
    plt.plot(x,y,'k-',linewidth=0.8*linew)
curve(2,7,0.8,0.69,-0.15),
text(2.5,3,'$y(x)$')
text(7.4,2.4,'$y(x+h)$')
text(6.5,3.9,'$y(x)+hf(x,y(x))$')
plt.show()
  • failure of exponential example
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .6
width = 12*scale
height = 8*scale
marker = 1.5*height
linew = marker/4
font = 3.5*height
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def line(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=0.8*linew)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/7,headwidth=marker,headlength=marker))
x0 = 1*scale
x1 = 12*scale
h = 0.1
A = -0.5
def curve(y0):
    x = np.linspace(x0,x1,100)
    y = scale*(4+(1+h*A)**(x/h)*y0)
    plt.plot(x,y,'-',color='0.75',linewidth=linew)
for y0 in range(-20,21,1):
    curve(y0)
def curvepoint(x,y0):
    x = x*scale
    y = scale*(4+(1+h*A)**(x/h)*y0)
    plt.plot(x,y,'ko',markersize=marker)
#for x in range(1,12,1):
#    curvepoint(x,4)
def step(h,N):
    x = x0
    h = h*scale
    y = 3
    for i in range(N):
        plt.plot(x,scale*(4+y),'ko',markersize=marker/1.2)
        xnew = x+h
        ynew = (1+h*A)*y
        plt.plot([x,xnew],[scale*(4+y),scale*(4+ynew)],'k--')
        x = xnew
        y = ynew
        plt.plot(x,scale*(4+y),'ko',markersize=marker/1.2)
step(1,10)
step(3,3)
step(5.4,2)
arrow(1,4,1,8)
arrow(1,4,1,0)
arrow(1,4,12,4)
text(12.25,4,'$x$')
text(.65,7,'$y$')
plt.show()
  • midpoint method
  • Runge-Kutta
  • step size error estimation

PDEs

  • 2D expansion
  • computational cluster
  • von Neumann stability
  • Courant-Friedrichs-Levy stability criterion
  • numerical dissipation
  • implicit methods
  • ADI
  • relaxation methods

Problems

  1. Solve the differential equation $\ddot x+x=0$ for an initial-value problem, and compare the dependence of the error at the end of the integration interval on the step size for an Euler and fourth-order Runge-Kutta method.
  2. Numerically solve for the motion of a parametrically driven pendulum $l\ddot\theta + (g + \ddot z)\sin \theta = 0$ and explore its dynamics.
  1. Numerically solved the wave equation $\pp{u}{t} = v^2\pp{u}{x}$, calculate its stability, and explore that boundary with and without artificial dissipation.