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¶
- 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.
- Numerically solve for the motion of a parametrically driven pendulum $l\ddot\theta + (g + \ddot z)\sin \theta = 0$ and explore its dynamics.
- 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.