Finite Differences
1. Airy equation¶
{\bf u}’ = \begin{bmatrix}0 & 1 \\ -t & 0 \end{bmatrix} {\bf u}
Copyright By PowCoder代写 加微信 powcoder
using LinearAlgebra, Plots, Interact
n = 50_000
t = range(0, 50; length=n)
A = t -> [0 1; -t 0]
h = step(t)
U = zeros(2,n) # each column is 𝐮_k
U[:,1] = [1,0] # initial condition
# Forward Euler:
for k = 1:n-1
U[:,k+1] = (I + h*A(t[k])) * U[:,k]
# for n = 50: we see unstable for Long time (though inaccurate)
plot(t, U’)
n = 100_000
t = range(0, 50; length=n)
A = t -> [0 1; -t 0]
h = step(t)
U = zeros(2,n) # each column is 𝐮_k
U[:,1] = [1,0] # initial condition
# Backward Euler:
for k = 1:n-1
U[:,k+1] = (I – h*A(t[k+1])) \ U[:,k] # solve each time step
# for n = 50: we see stable for Long time (though inaccurate)
plot(t, U’)
2. Heat on a Graph¶
{\bf u}’ = \underbrace{\begin{bmatrix}-1 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & 1 \\ && 1 & -2 & 1 \\ &&&1 & -1 \end{bmatrix}}_Δ {\bf u} + e_{floor(m/2)} * cos(ω*t)
m = 50 # 50 nodes
SymTridiagonal([-1.; fill(-2, m-2); -1.], fill(1., m-1))
U = zeros(m,n) # each column is 𝐮_k
U[:,1] = zeros(m) # initial condition
𝐞 = zeros(m); 𝐞[m÷2] = 1;
# Forward Euler:
for k = 1:n-1
U[:,k+1] = (I + h*Δ) * U[:,k] + h * 𝐞 * cos(ω*t[k])
# Backward Euler:
for k = 1:n-1
U[:,k+1] = (I – h*Δ) \ (U[:,k] + h * 𝐞 * (2 + cos(ω*t[k+1])))
@manipulate for k=1:n
scatter(U[:,k]; ylims=(-1,5), label=”time = $((k-1)*h)”)
3. Pendulum¶
{\bf u}’ = \begin{bmatrix}u_2(t) \\ -\sin(u_1(t)) \end{bmatrix}
using LinearAlgebra, Plots, Interact
n = 50_000
t = range(0, 50; length=n)
h = step(t)
U = zeros(2,n) # each column is 𝐮_k
U[:,1] = [0,1.6] # initial condition
# Forward Euler:
for k = 1:n-1
U[:,k+1] = U[:,k] + h * [U[2,k], -sin(U[1,k])]
plot(t, U’)
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com