In [1]:
using Plots,LaTeXStrings
default(markersize=3,linewidth=1.5)
using LinearAlgebra,DifferentialEquations
using SparseArrays
include("FNC.jl");
Consider again the example $f(x,y)=\sin(xy)$ over $[0,2\pi]\times[1,3]$, with the dimensions discretized using $m=4$ and $n=2$ equal pieces, respectively.
In [2]:
m = 4; x = 0:2*pi/m:2*pi;
n = 2; y = 1:2/n:3;
We create a representation of the grid using two matrices created by the ndgrid function.
In [3]:
X,Y = FNC.ndgrid(x,y);
X
Out[3]:
In [4]:
Y
Out[4]:
As you see above, the entries of X vary in the first dimension (rows), while the entries of Y vary along the second dimension (columns).
We can also visualize this grid on the rectangle.
In [5]:
scatter(X,Y,m=:blue,xtick=x,ytick=y,leg=:none)
Out[5]:
For a given definition of $f(x,y)$ we can find $\operatorname{mtx}(f)$ by elementwise operations on the coordinate matrices X and Y.
In [6]:
f = (x,y) -> sin(x*y-y);
F = f.(X,Y)
Out[6]:
We can make nice plots of the function by first choosing a much finer grid.
In [7]:
m = 70; x = 0:2*pi/m:2*pi;
n = 50; y = 1:2/n:3;
X,Y = FNC.ndgrid(x,y);
F = f.(X,Y);
plot(x,y,F,levels=10,xlabel=L"x",ylabel=L"y",fill=true,match_dimensions=true)
Out[7]:
In [8]:
surface(x,y,F,xlabel="x",ylabel="y",zlabel="f(x,y)",match_dimensions=true,leg=:false)
Out[8]:
Construction of a function over the unit disk is straightforward. We define a grid in $(r,\theta)$ space and compute accordingly. For the function $f(r,\theta)=1-r^4$, for example:
In [9]:
r = (0:40)/40
theta = 2pi*(0:80)/80
R,Theta = FNC.ndgrid(r,theta);
F = @. 1-R^4;
In [10]:
surface(r,theta,F,match_dimensions=true,
xlabel="r",ylabel="\\theta",title="A polar function",leg=:none)
Out[10]:
Of course we are used to seeing such plots over the $(x,y)$ plane. To do this in Plots at this writing, we need to switch the graphics backend. (These graphs do not appear in the plain HTML version of the page. To see them you must view the notebook source.)
In [11]:
plotlyjs()
X = @. R*cos(Theta); Y = @. R*sin(Theta);
surface(X,Y,F,xlabel="x",ylabel="y",leg=:none)
Out[11]:
In such functions it's up to us to ensure that the values along the line $r=0$ are identical, and that the values on the line $\theta=0$ are identical to those on $\theta=2\pi$. Otherwise the interpretation of the domain as the unit disk is nonsensical.
On the unit sphere we can use color to indicate a function value. For the function $f(x,y,z)=xyz^3$, for instance:
In [12]:
theta = 2pi*(0:60)/60
phi = pi*(0:60)/60
Theta,Phi = FNC.ndgrid(theta,phi);
X = @.cos(Theta)*sin(Phi)
Y = @.sin(Theta)*sin(Phi)
Z = cos.(Phi)
F = @. X*Y*Z^3
surface(X,Y,Z,fill_z=F,
xlabel="x",ylabel="y",zlabel="z",title="Function on the unit sphere")
Out[12]:
In [13]:
gr();
We define a function and, for reference, its two exact partial derivatives.
In [14]:
u = (x,y) -> sin(x*y-y);
dudx = (x,y) -> y*cos(x*y-y);
dudy = (x,y) -> (x-1)*cos(x*y-y);
We use an equispaced grid and second order finite differences as implemented by diffmat2.
In [15]:
m = 80; x,Dx = FNC.diffmat2(m,[0,2*pi]);
n = 60; y,Dy = FNC.diffmat2(n,[1,3]);
X,Y = FNC.ndgrid(x,y)
U = u.(X,Y);
Here we compare the exact $\partial u/\partial x$ with the finite-difference approximation.
In [16]:
contourf(x,y,dudx.(X,Y),match_dimensions=true,levels=10,layout=(1,2),subplot=1,
title="du/dx")
contourf!(x,y,Dx*U,levels=10,match_dimensions=true,layout=(1,2),subplot=2,
title="approximation")
Out[16]:
We now do the same for $\partial u/\partial y$,
In [17]:
contourf(x,y,dudy.(X,Y),match_dimensions=true,levels=10,layout=(1,2),subplot=1,
title="du/dy")
contourf!(x,y,U*Dy',levels=10,match_dimensions=true,layout=(1,2),subplot=2,
title="approximation")
Out[17]:
To the eye there is little difference to be seen, though we expect that the results probably have no more than 2-4 correct digits for these discretization sizes.
In [18]:
m = 2; n = 3;
v = [1:6;]
Out[18]:
The unvec operation:
In [19]:
V = reshape(v,2,3)
Out[19]:
Two equivalent syntaxes for the vec operation:
In [20]:
vec(V)
Out[20]:
In [21]:
V[:]
Out[21]:
We will solve a 2D heat equation $u_t = 0.1(u_{xx} + u_{yy})$ on the square $[-1,1]\times[-1,1]$. We'll assume periodic behavior in both directions.
In [22]:
m = 60; x,Dx,Dxx = FNC.diffper(m,[-1,1]);
n = 40; y,Dy,Dyy = FNC.diffper(n,[-1,1]);
X,Y = FNC.ndgrid(x,y);
Note that the initial condition must also be periodic on the domain.
In [23]:
U0 = @. sin(4*pi*X)*exp(cos(pi*Y))
contourf(x,y,U0,match_dimensions=true,color=:redsblues,aspect_ratio=1,
xaxis=(L"x"),yaxis=(L"y"),title="Initial condition")
Out[23]:
The next two functions map between the natural matrix shape of the unknowns and the vector shape demanded by the ODE solvers.
In [24]:
function unpack(u) reshape(u,m,n); end
function pack(U) U[:]; end
Out[24]:
This function computes the time derivative for the unknowns. The actual calculations take place using the matrix shape.
In [25]:
function dudt(u,nu,t)
U = unpack(u);
Uxx = Dxx*U; Uyy = U*Dyy'; # 2nd partials
dUdt = nu*(Uxx + Uyy); # PDE
return pack(dUdt);
end
Out[25]:
Since this problem is parabolic, a stiff integrator is a good choice.
In [26]:
IVP = ODEProblem(dudt,pack(U0),(0,0.2),0.1)
sol = solve(IVP,Rodas4P());
Here we plot the solution at two different times. (The results are best viewed using an animation.)
In [27]:
an = @animate for t = range(0,stop=0.2,length=81)
surface(x,y,unpack(sol(t)),match_dimensions=true,color=:redsblues,clims=(-2,2),
xaxis=(L"x"),yaxis=(L"y"),title="Heat equation, t=$(round(t,digits=3))")
end
gif(an,"heat2Dper.gif")
Out[27]:
We will solve an advection-diffusion problem, $u_t + u_x = 1 + \nu(u_{xx} + u_{yy})$, where $u=0$ on the boundary of the square $[-1,1]^2$.
In [28]:
m = 50; n = 50;
X,Y,d = FNC.rectdisc(m,[-1,1],n,[-1,1]);
The initial condition we specify here is used to impose its boundary values on the solution at all times.
In [29]:
U0 = @. (1-X^4)*(1-Y^4);
plot(X[:,1],Y[1,:],U0,match_dimensions=true,
xlabel=L"x",ylabel=L"y",aspect_ratio=1,title="Initial condition")
Out[29]:
This next function maps the unknowns, given in a vector shape, to a matrix of values including the boundaries.
In [30]:
function unpack(w)
U = copy(U0) # get the boundary right
U[@. !d.isbndy] .= w # overwrite the interior
return U
end
Out[30]:
The next function drops the boundary values and returns a vector of the interior values. It's the inverse of the unpack function.
In [31]:
function pack(U) U[@. !d.isbndy]; end
Out[31]:
This function computes the time derivative at the interior nodes only.
In [32]:
function dwdt(w,nu,t)
U = unpack(w)
Uxx = d.Dxx*U; Uyy = U*d.Dyy'; # 2nd partials
dUdt = 1 .- d.Dx*U + nu*(Uxx + Uyy); # PDE
return pack(dUdt)
end
Out[32]:
Since this problem is parabolic, a stiff integrator is a good choice. (The solver usually does fine figuring this out on its own, but here we give it a hint.)
In [33]:
IVP = ODEProblem(dwdt,pack(U0),(0.0,1.0),0.05)
sol = solve(IVP,alg_hints=[:stiff]);
In [34]:
plot(X[:,1],Y[1,:],unpack(sol(0.5)),match_dimensions=true,
xlabel=L"x",ylabel=L"y",aspect_ratio=1,title="Solution at t=0.5")
Out[34]:
In [35]:
plot(X[:,1],Y[1,:],unpack(sol(1)),match_dimensions=true,
xlabel=L"x",ylabel=L"y",aspect_ratio=1,title="Solution at t=1")
Out[35]:
In [36]:
an = @animate for t in range(0,stop=1,length=80)
surface(X[:,1],Y[1,:],unpack(sol(t)),match_dimensions=true,aspect_ratio=1,color=:blues,clims=(0,2),
title="Advection-diffusion solution at t=$(round(t,digits=2))")
end
gif(an,"advdiff2D.gif")
Out[36]:
In [37]:
m = 60; n = 60;
X,Y,d = FNC.rectdisc(m,[-2,2],n,[-2,2]);
x = X[:,1]; y = Y[1,:]; # for plotting
Here is the initial condition. The boundary values of $u$ will remain constant.
In [38]:
U0 = @. (X+0.2)*exp(-12*(X^2+Y^2))
V0 = zeros(size(U0))
surface(x,y,U0,
title="Initial condition",xaxis=(L"x"),yaxis=(L"y") )
Out[38]:
The unpack function separates the unknowns for $u$ and $v$, applies the boundary conditions on $u$, and returns two functions on the grid.
In [39]:
function unpack(w)
numU = (m-1)*(n-1) # number of unknowns for U
U = copy(U0)
U[@. !d.isbndy] = w[1:numU] # overwrite the interior
V = d.unvec( w[numU+1:end] ) # use all values
return U,V
end
Out[39]:
The next function drops the boundary values of $u$ and returns a vector of all the unknowns for both components of the solution. It's the inverse of the unpack function.
In [40]:
function pack(U,V)
w = U[@. !d.isbndy]
return [ w; V[:] ]
end
Out[40]:
The following function computes the time derivative of the unknowns. Besides the translation between vector and matrix shapes, it's quite straightforward.
In [41]:
function dwdt(w,tmp,t)
U,V = unpack(w)
dUdt = V
dWdt = d.Dxx*U + U*d.Dyy'
return pack(dUdt,dWdt)
end
Out[41]:
Since this problem is hyperbolic, not parabolic, a nonstiff integrator like |ode45| is fine and faster than a stiff integrator.
In [42]:
IVP = ODEProblem(dwdt,pack(U0,V0),(0,3.0))
sol = solve(IVP,alg_hints=[:nonstiff]);
In [43]:
an = @animate for t = range(0,stop=3,length=80)
surface(x,y,unpack(sol(t)),match_dimensions=true,color=:redsblues,clims=(-0.1,0.1),
title="Wave equation solution at t=$(round(t,digits=2))")
end
gif(an,"wave2D.gif")
Out[43]:
In [44]:
A = [1 2; -2 0]
B = [ 1 10 100; -5 5 3 ]
Out[44]:
Applying the definition manually,
In [45]:
A_kron_B = [ A[1,1]*B A[1,2]*B;
A[2,1]*B A[2,2]*B ]
Out[45]:
But it's simpler to use the built-in kron.
In [46]:
kron(A,B)
Out[46]:
Here is a forcing function for Poisson's equation.
In [47]:
phi = (x,y) -> x^2 - y + 2;
We pick a crude discretization for illustrative purposes.
In [48]:
m = 5; n = 6;
X,Y,d = FNC.rectdisc(m,[0,3],n,[-1,1]);
Next, we evaluate $\phi$ on the grid.
In [49]:
F = phi.(X,Y);
Here are the equations for the PDE collocation, before any modifications are made for the boundary conditions.
In [50]:
A = kron(d.Iy,d.Dxx) + kron(d.Dyy,d.Ix);
spy(sparse(A),color=:reds,title="Matrix before BC")
Out[50]:
In [51]:
b = d.vec(F);
@show N = length(b);
The array d.isbndy is Boolean and the same size as X, Y, and F.
In [52]:
spy(sparse(d.isbndy),m=5,color=:reds,title="Boundary points")
plot!([],label="",grid=:xy,xaxis=([0,8]),yaxis=([0,7]))
Out[52]:
Next replace the boundary rows of the system by rows of the identity.
In [53]:
I = spdiagm(0=>ones(size(A,1)))
A[d.isbndy[:],:] = I[d.isbndy[:],:] # Dirichlet conditions
spy(sparse(A),color=:reds,title="Matrix with BC")
Out[53]:
Finally, we must replace the rows in the vector $\mathbf{b}$ by the boundary values being assigned to the boundary points. Here, we let the boundary values be zero everywhere.
In [54]:
b[d.isbndy[:]] .= 0; # Dirichlet values
Now we can solve for $\mathbf{u}$ and reinterpret it as the matrix-shaped $\mathbf{U}$, the solution on our grid. This grid is much too coarse for the result to look like a smooth function of two variables.
In [55]:
u = A\b;
U = d.unvec(u);
wireframe(X[:,1],Y[1,:],U,match_dimensions=true,
xaxis=("x"),yaxis=("y"),zaxis=("u(x,y)"),title="Coarse solution")
Out[55]:
First we define the problem on $[0,1]\times[0,2]$.
In [56]:
f = (x,y) -> -sin(3*x.*y-4*y)*(9*y^2+(3*x-4)^2);
g = (x,y) -> sin(3*x*y-4*y);
xspan = [0,1]; yspan = [0,2];
Here is the finite difference solution.
In [57]:
U,X,Y = FNC.poissonfd(f,g,50,xspan,80,yspan);
In [58]:
x = X[:,1]; y = Y[1,:];
surface(x,y,U,match_dimensions=true,color=:bluesreds,
title="Solution of Poisson equation",
xaxis=("x"), yaxis=("y"), zaxis=("u(x,y)"))
Out[58]:
The error is a surprisingly smooth function.
In [59]:
error = g.(X,Y) - U;
M = max( maximum(error),-minimum(error) )
contourf(x,y,error,levels=17,match_dimensions=true,aspect_ratio=1,clims=(-M,M),color=:bluesreds,
title="Error",xaxis=("x"),yaxis=("y") )
Out[59]:
Because $u$ is specified on the boundary, we should observe that the error is zero on the boundary---else we have done something incorrectly.
The function here defines the PDE, by way of $\mathbf{f}(\mathbf{u})$ and its derivative (Jacobian matrix). Note that the function will have access to all of the properties of a discretization, as if they were returned by rectdisc.
In [60]:
lambda = 1.5
function pde(U,X,Y,d)
LU = d.Dxx*U + U*d.Dyy'; # apply Laplacian
F = @. LU - lambda/(U+1)^2 # residual
L = kron(d.Dyy,d.Ix) + kron(d.Iy,d.Dxx)
u = d.vec(U)
J = L + spdiagm( 0 => @. 2*lambda/(u+1)^3 )
return F,J
end
Out[60]:
We define a trivial zero function for the Dirichlet boundary condition, and solve.
In [61]:
g = (x,y) -> 0 # boundary condition
U,X,Y = FNC.newtonpde(pde,g,100,[0,2.5],80,[0,1]);
In [62]:
surface(X[:,1],Y[1,:],U',color=:blues,aspect_ratio=1,
xlabel="x",ylabel="y",title="Deflection of a MEMS membrane")
Out[62]:
The following defines the PDE, by way of a function for $\mathbf{f}(\mathbf{u})$ and its derivative (Jacobian matrix), and a function for the Dirichlet boundary condition.
In [63]:
function pde(U,X,Y,d)
LU = d.Dxx*U + U*d.Dyy' # apply Laplacian
F = @. U*(1-U^2) + 0.05*LU # residual
L = kron(d.Dyy,d.Ix) + kron(d.Iy,d.Dxx)
u = d.vec(U)
J = spdiagm(0 => @. 1-3*u^2) + 0.05*L # Jacobian
return F,J
end
g = (x,y) -> tanh(5*(x+2*y-1)); # boundary condition
We solve the PDE and then plot the result.
In [64]:
U,X,Y = FNC.newtonpde(pde,g,100,[0,1],100,[0,1]);
In [65]:
surface(X[:,1],Y[1,:],U',camera=(32,30),xlabel="x",ylabel="y",
title="Steady Allen-Cahn")
Out[65]:
The following defines the PDE, by way of a function for $\mathbf{f}(\mathbf{u})$ and its derivative (Jacobian matrix), and a function for the Dirichlet boundary condition.
In [66]:
function pde(U,X,Y,d)
LU = d.Dxx*U + U*d.Dyy' # apply Laplacian
Ux = d.Dx*U
F = @. 1 - Ux + 0.05*LU # residual
L = kron(d.Dyy,d.Ix) + kron(d.Iy,d.Dxx)
u = d.vec(U)
J = -kron(d.Iy,d.Dx) + 0.05*L # Jacobian
return F,J
end
g = (x,y) -> 0; # boundary condition
In [67]:
U,X,Y = FNC.newtonpde(pde,g,100,[-1,1],100,[-1,1]);
In [68]:
surface(X[:,1],Y[1,:],U',
xlabel="x",ylabel="y",title="Steady Advection-diffusion")
Out[68]: