{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modules and packages\n",
"### Default library\n",
"Julia has a large number of default functionalities that can be used for many purposes in generic and scientific computations. Modules in default library can be loaded with `using ModuleName`. Let's check out some of the functionalities of the `LinearAlgebra` module, for example:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`LinearAlgebra` module provides native implementations of many common and useful linear algebra operations. For example, the inverse of a matrix `A` can be easily computed with `inv(A)`:"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.5 -2.0 1.5\n",
" -2.0 1.0 0.0\n",
" 1.5 0.0 -0.5"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 2 3;2 5 6;3 6 7]\n",
"inv(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or its eigenvalues can be easily computed with `eigvals(A)`:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" -0.46672160014668007\n",
" 0.32610414765071355\n",
" 13.140617452495963"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Complex operations/factorizations are of course supported (e.g., skew-symmetric matrix, 90° planar rotation):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}\n",
"values:\n",
"2-element Array{Complex{Float64},1}:\n",
" 0.0 - 1.0im\n",
" 0.0 + 1.0im\n",
"vectors:\n",
"2×2 Array{Complex{Float64},2}:\n",
" 0.707107-0.0im 0.707107+0.0im\n",
" 0.0+0.707107im 0.0-0.707107im"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigen([0 -1;1 0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`LinearAlgebra` does not only define new functions, it defines new types for the most common matrix structures, such as `Symmetric`, `Hermitian`, `UpperTriangular`, `LowerTriangular`, `Diagonal`, `Tridiagonal`, `SymTridiagonal`, `Bidiagonal`, `UniformScaling`. These type definitions, together with methods definitions, follow the philosophy of multiple dispatch: let's define a symmetric matrix `A` and let's create a `Symmetric` instantiation `symm_A` from it:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Symmetric{Int64,Array{Int64,2}}:\n",
" -2 1\n",
" 1 -2"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [-2 1;1 -2]\n",
"symm_A = Symmetric(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`A` is treated as a generic matrix since we didn't create it with a particular type, therefore any matrix operation (e.g., multiplication) on it will employ a generic method for the supertype `AbstractArray` for generic arrays:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"*(A::AbstractArray{T,2} where T , B::AbstractArray{T,2} where T ) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\matmul.jl:151 "
],
"text/plain": [
"*(A::AbstractArray{T,2} where T, B::AbstractArray{T,2} where T) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\matmul.jl:151"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@which A*A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"while the same matrix multiplication with a `Symmetric` type will use a different and more specific method:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"*(A::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T , B::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T ) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\symmetric.jl:569 "
],
"text/plain": [
"*(A::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T, B::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\symmetric.jl:569"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@which symm_A*symm_A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another quick example is the `tr()` operator $\\text{tr}(\\mathbf{A})=\\sum a_{ii}$, with different methods for different types of the input matrix:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"# 7 methods for generic function tr :"
],
"text/plain": [
"# 7 methods for generic function \"tr\":\n",
"[1] tr(A::Array{T,2}) where T in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\dense.jl:330\n",
"[2] tr(A::Hermitian) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\symmetric.jl:363\n",
"[3] tr(D::Diagonal) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\diagonal.jl:562\n",
"[4] tr(A::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti} where Ti<:Integer) where Tv in SparseArrays at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\SparseArrays\\src\\sparsematrix.jl:3530\n",
"[5] tr(A::AbstractArray{T,2} where T) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\generic.jl:1007\n",
"[6] tr(x::Number) in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\generic.jl:1011\n",
"[7] tr(J::UniformScaling{T}) where T in LinearAlgebra at C:\\Users\\riki\\AppData\\Local\\JuliaPro-1.5.3-1\\Julia-1.5.3\\share\\julia\\stdlib\\v1.5\\LinearAlgebra\\src\\uniformscaling.jl:215"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"methods(tr)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This implementation of `LinearAlgebra` (as all other modules) follows the multiple dispatch philosophy where different methods are employed for the same function, depending upon the argument types. In this way the code can be easily extended: suppose we want to introduce new functionalities for penta-diagonal matrices. All we need to do is to define a new type `Pentadiagonal`, which will be a subtype of `AbstractArray`, and add new method definitions to perform the required matrix operations on this new array type only, using multiple dispatch. Existing type/methods definitions do not need any modification.\n",
"\n",
"Usual matrix factorizations such as `Cholesky`, `LDLt`, `LU`, `QR`, `SVD`, etc. are available:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"LU{Float64,Array{Float64,2}}\n",
"L factor:\n",
"2×2 Array{Float64,2}:\n",
" 1.0 0.0\n",
" -0.5 1.0\n",
"U factor:\n",
"2×2 Array{Float64,2}:\n",
" -2.0 1.0\n",
" 0.0 -1.5"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ch = lu(symm_A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`factorize(A)` compute a convenient factorization of A, based upon the type of the input matrix. For example, in the case of a symmetric (but not definite positive) matrix, a symmetric indefinite `BunchKaufman` factorization ($\\mathbf{A}=\\mathbf{L}\\mathbf{D}\\mathbf{L}^T$ or $\\mathbf{A}=\\mathbf{U}\\mathbf{D}\\mathbf{U}^T$ with permutation) is automatically chosen:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BunchKaufman{Float64,Array{Float64,2}}\n",
"D factor:\n",
"2×2 Tridiagonal{Float64,Array{Float64,1}}:\n",
" -1.5 0.0\n",
" 0.0 -2.0\n",
"U factor:\n",
"2×2 UnitUpperTriangular{Float64,Array{Float64,2}}:\n",
" 1.0 -0.5\n",
" ⋅ 1.0\n",
"permutation:\n",
"2-element Array{Int64,1}:\n",
" 1\n",
" 2"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"factorize(symm_A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"while in the case of a non symmetric `B`, a `LU` factorization ($\\mathbf{A}=\\mathbf{L}\\mathbf{U}$) is returned:"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"LU{Float64,Tridiagonal{Float64,Array{Float64,1}}}\n",
"L factor:\n",
"2×2 Array{Float64,2}:\n",
" 1.0 0.0\n",
" 0.333333 1.0\n",
"U factor:\n",
"2×2 Array{Float64,2}:\n",
" 3.0 4.0\n",
" 0.0 0.666667"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = [1 2;3 4]\n",
"LU_B = factorize(B)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Factorizations can be directly employed when the solution of a linear systems is required, as if it was the original matrix:"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.9999999999999998\n",
" -0.49999999999999983"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"RHS = [1;4]\n",
"x = B\\RHS"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.9999999999999998\n",
" -0.49999999999999983"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#B * x = RHS\n",
"x = LU_B\\RHS"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where in the last case the LU factorization is actually employed, saving computations for multiple RHS."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Conjugate transpose (adjoint) is performed with the usual `'`. Note that the resulting array object has its own type `Adjoint`:"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Adjoint{Int64,Array{Int64,2}}:\n",
" 1 3\n",
" 2 4"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 2;3 4]\n",
"B = A'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where `B` is a \"lazy\" adjoint, i.e., it is not a fully independent array storing all the elements of the adjoint of `A`. Therefore if the entries of `A` are modified, the corresponding entries of `B` are also modified:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Adjoint{Int64,Array{Int64,2}}:\n",
" 100 3\n",
" 2 4"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A[1] = 100\n",
"B"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`Adjoint` is again a (parametric) composite type, therefore we can check where `B::Adjoint` is pointing by checking the name of the field(s) of this composite type:"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(:parent,)"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fieldnames(Adjoint)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`parent` is therefore the only field that is needed to lazily perform an adjoint. We can check if the `parent` of `B` points to `A`:"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ptr{Nothing} @0x000000003a2972d0\n",
"Ptr{Nothing} @0x000000003a2972d0\n"
]
}
],
"source": [
"println( pointer_from_objref(A) )\n",
"println( pointer_from_objref(B.parent) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Packages\n",
"Additional functionalitites which are not present in the default library can be added to Julia by means of packages. Any external package `PackageName` can be added by typing:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Pkg\n",
"Pkg.add(\"PackageName\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Plots\n",
"`Plots` (\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = range(0, 2pi, length = 100)\n",
"p1 = plot(x, sin)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Multiple plots can be displayed without being overwritten (`hold on` in MATLAB) by using the `!` name convention for functions that modify their arguments. In this case the new plot is added to `p1`:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"display( plot!(p1, x, cos) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plots are automatically displayed in interactive sessions (i.e., in the REPL or here in Jupyter), while `display(plot(...))` is needed to display plots in scripts.\n",
"\n",
"By using attributes defined by keywords, we can customize plot appearance (more on attributes here \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0 \n",
" \n",
"\n",
"1 \n",
" \n",
"\n",
"2 \n",
" \n",
"\n",
"3 \n",
" \n",
"\n",
"4 \n",
" \n",
"\n",
"5 \n",
" \n",
"\n",
"6 \n",
" \n",
"\n",
"−1.0 \n",
" \n",
"\n",
"−0.5 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"1.0 \n",
" \n",
" \n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x, sin, w = 3, c = :red,\n",
" tickfont = (:times,10),\n",
" ylabel = \"\\$y\\$\", xlim = (0,2pi),\n",
" xlabel = \"\\$x\\$\", ylim = (-1,1),\n",
" guidefontsize = 14,\n",
" size = (400,200),\n",
" leg = false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Figures can be exported to standard formats such as `pdf`, `png`, `eps`, `ps`, `svg`, `html`, `tex`, depending on the backend:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [],
"source": [
"savefig(\"myplot.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Different types of plots are supported (2D, 3D, etc.). Let's check out `heatmap` which can be employed to display bivariate functions $f(x,y)$ as color maps:"
]
},
{
"cell_type": "code",
"execution_count": 113,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
" \n",
" \n"
]
},
"execution_count": 113,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = range(0, 1, length = 500)\n",
"y = range(0, 1, length = 500)\n",
"f(x,y) = sin(2pi/(x+.1)) * sin(2pi*y)\n",
"heatmap(x, y, f,\n",
" xlim = (0,1), ylim = (0, 1), size=(900,400),\n",
" colorbar = false,\n",
" tickfont = (:times,12))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Animations can be done very easily by using the macro `@gif`, for example:"
]
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = C:\\Users\\riki\\tmp.gif\n",
"└ @ Plots C:\\Users\\riki\\.julia\\packages\\Plots\\SVksJ\\src\\animation.jl:104\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
"\n",
"0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.7 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"0.9 \n",
" \n",
"\n",
"1.0 \n",
" \n",
" \n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 100\n",
"φ = zeros(N, N)\n",
"i = j = 2:(N-1)\n",
"φ[i,j] .= rand(N-2, N-2)\n",
"\n",
"Δx = 1 / (N-1)\n",
"Δt = Δx^2 / 5\n",
"for n ∈ 0:500\n",
" φ[i,j] += Δt * ( φ[i.+1,j] + φ[i.-1,j] + φ[i,j.+1] + φ[i,j.-1] - 4 * φ[i,j] ) / Δx^2\n",
"end\n",
"\n",
"x = y = range(0, 1, length = N)\n",
"heatmap(x, y, φ, xlim = (0,1), ylim = (0, 1), aspect_ratio = 1, tickfont = (:times,10), clim=(0,1)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If $\\phi$ is interpreted as a temperature field, we are solving an unsteady heat conduction problem on a square domain whose edges are constantly kept \"cool\" at $\\phi=0$. We can enforce a hot spot inside $\\Omega$ (e.g., constant \"hot\" temperature $\\phi=1$) and create an animation to figure out which is the time-dependent temperature distribution during the heating:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
"\n",
"0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.7 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"0.9 \n",
" \n",
"\n",
"1.0 \n",
" \n",
" \n"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function explicit_for_loop_SOR!(A; n_iterations, ω)\n",
" N = size(A, 1)\n",
" k = 2:(N-1)\n",
" hot_spot_location = (125, 100:400)\n",
" for iteration ∈ 1:n_iterations\n",
" φ[hot_spot_location...] .= 1\n",
" for i ∈ k, j ∈ k\n",
" A[i,j] += ω * ( ( A[i+1,j] + A[i-1,j] + A[i,j+1] + A[i,j-1] ) / 4 - A[i,j] )\n",
" end\n",
" end\n",
" return nothing\n",
"end\n",
"\n",
"φ .= 0\n",
"@time explicit_for_loop_SOR!(φ, n_iterations=1000, ω=1.5)\n",
"\n",
"x = y = range(0, 1, length = N)\n",
"h1 = heatmap(x, y, φ,\n",
" xlabel = \"\\$x\\$\", xlim = (0,1),\n",
" ylabel = \"\\$y\\$\", ylim = (0,1),\n",
" tickfont = (:times,10), aspect_ratio = 1, clim=(0,1)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With `Plots` we can display very easily some temperature profiles along $y$. We can also use module `Printf` for formatted output:"
]
},
{
"cell_type": "code",
"execution_count": 113,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.00 \n",
" \n",
"\n",
"0.25 \n",
" \n",
"\n",
"0.50 \n",
" \n",
"\n",
"0.75 \n",
" \n",
"\n",
"1.00 \n",
" \n",
"\n",
"0.00 \n",
" \n",
"\n",
"0.25 \n",
" \n",
"\n",
"0.50 \n",
" \n",
"\n",
"0.75 \n",
" \n",
"\n",
"1.00 \n",
" \n",
"\n",
" \n"
]
},
"execution_count": 113,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Printf\n",
"v = 1:51:N\n",
"φ_profiles = φ[:,v]\n",
"plot(y, φ_profiles,\n",
" label = [ @sprintf(\"\\$x=%.2f\\$\", x[j]) for i ∈ [1], j ∈ v ],\n",
" legendfontsize = 10,\n",
" tickfont = (:times,10),\n",
" xlabel = \"\\$y\\$\",\n",
" ylabel = \"\\$T\\$\",\n",
" palette = palette(:thermal, length(v)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also evaluate numerically the gradient of our computed field `φ` in a 2-component (named) tuple of 2D arrays `∇φ`:\n",
"\n",
"$\\qquad\\displaystyle\\frac{\\partial \\phi}{\\partial x}\\Big\\vert_{i,j}\\approx\\frac{\\phi_{i+1,j}-\\phi_{i-1,j}}{2\\Delta x}$\n",
"\n",
"$\\qquad\\displaystyle\\frac{\\partial \\phi}{\\partial y}\\Big\\vert_{i,j}\\approx\\frac{\\phi_{i,j+1}-\\phi_{i,j-1}}{2\\Delta y}$\n",
"\n",
"in order to put the vector plot of the flux on the previous heatmap `h1`:"
]
},
{
"cell_type": "code",
"execution_count": 188,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
" \n",
"\n",
" \n",
"\n",
"0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.7 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"0.9 \n",
" \n",
"\n",
"1.0 \n",
" \n",
" \n"
]
},
"execution_count": 188,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"∇φ = ( ∂x=(φ[i.+1,j] .- φ[i.-1,j])/(2*Δx) , ∂y=(φ[i,j.+1] .- φ[i,j.-1])/(2*Δx) )\n",
"v = 19:19:(N-2)\n",
"x_in = x[i][v]\n",
"X = (1 .+ 0*x_in) * x_in'\n",
"𝜆 = 0.01\n",
"quiver!(h1, X[:], X'[:],\n",
" quiver = ( -𝜆*∇φ.∂y[v,v][:] , -𝜆*∇φ.∂x[v,v][:] ),\n",
" xlim = (0,0.5), ylim = (0,0.5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Parallelization\n",
"Can we do even better in terms of performance by using explicit `for` loops? If our machine has a multi-core processor (which is extremely likely), we can distribute our computations so that different calculations are executed simultaneously on each core. In order to exploit parallelization on such shared-memory architectures, the environment variable `JULIA_NUM_THREADS` must be set to the number `n` of cores/physical threads of your CPU (e.g., `n=2, 4, 8, 16`) before starting Julia:\n",
"\n",
"* Windows: in command prompt, type `setx JULIA_NUM_THREADS n`\n",
"* Unix: in bash, type `export JULIA_NUM_THREADS=n`\n",
"\n",
"If you don't know the number of cores `n` of your CPU, in Windows, type `echo %NUMBER_OF_PROCESSORS%` in command prompt.\n",
"\n",
"Then, we need to load the required module for parallelization, which is already present in the `Base` library:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"using Base.Threads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If `JULIA_NUM_THREADS` has been set correctly, now we should get the number `n`:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nthreads()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the `@threads` macro to equally subdivide the iterations of a `for` loop over the initialized threads:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 100\n",
"v = zeros(Int64, N)\n",
"#v = (1:N)/N\n",
"@threads for i ∈ eachindex(v)\n",
" v[i] = threadid()\n",
" # some other computation\n",
"end\n",
"scatter(v, yticks = 1:8, legend = false)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 5.624906 seconds (555.30 k allocations: 73.311 MiB, 0.35% gc time)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
" \n",
"\n",
" \n",
"\n",
"0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.7 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"0.9 \n",
" \n",
"\n",
"1.0 \n",
" \n",
" \n"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function explicit_parallel_for_loop_SOR!(A; n_iterations, ω)\n",
" N = size(A, 1)\n",
" k = 2:(N-1)\n",
" hot_spot_location = (125, 100:400)\n",
" for iteration ∈ 1:n_iterations\n",
" φ[hot_spot_location...] .= 1\n",
" @threads for i ∈ k\n",
" for j ∈ k\n",
" A[i,j] += ω * ( ( A[i+1,j] + A[i-1,j] + A[i,j+1] + A[i,j-1] ) / 4 - A[i,j] )\n",
" end\n",
" end\n",
" end\n",
" return nothing\n",
"end\n",
"\n",
"φ .= 0\n",
"@time explicit_parallel_for_loop_SOR!(φ, n_iterations=10000, ω=1.5)\n",
"\n",
"x = y = range(0, 1, length = N)\n",
"h2 = heatmap(x, y, φ,\n",
" xlabel = \"\\$x\\$\", xlim = (0,1),\n",
" ylabel = \"\\$y\\$\", ylim = (0,1),\n",
" tickfont = (:times,10), aspect_ratio = 1, clim=(0,1)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Solution via direct method\n",
"Instead of computing an approximation of the steady-state solution through an iterative method, in Julia we can easily compute the actual steady solution by solving the corresponding linear system:\n",
"\n",
"$4\\phi_{i,j}-\\phi_{i+1,j}-\\phi_{i-1,j}-\\phi_{i,j+1}-\\phi_{i,j-1}=0\\qquad\\big(i.e.,4\\phi_{node}-\\sum\\phi_{neighbour\\;nodes}=0\\big)$\n",
"\n",
"for each of the $M^2$ internal nodes. Let's create the coefficients $a$ for this linear system, created as $M\\times M$ arrays:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"M = 500\n",
"a_P = 4*ones(M, M) # i , j\n",
"a_E = -ones(M, M) # i+1 , j\n",
"a_W = -ones(M, M) # i-1 , j\n",
"a_N = -ones(M, M) # i , j+1\n",
"a_S = -ones(M, M) # i , j-1\n",
"\n",
"# BCs\n",
"a_E[M,:] .= 0\n",
"a_W[1,:] .= 0\n",
"a_N[:,M] .= 0\n",
"a_S[:,1] .= 0 ;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have to actually build the $M^2\\times M^2$ matrix $\\mathbf{A}$ of the corresponding linear system. $\\mathbf{A}$ is sparse since each row has only 5 nonzero coefficients, therefore we can use the `SparseArrays` module, already present in the standard library, which allows easy creation/manipulation of sparse matrices. `spdiagm` creates a sparse matrix from its nonzero diagonals:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"250000×250000 SparseMatrixCSC{Float64,Int64} with 1248998 stored entries:\n",
" [1 , 1] = 4.0\n",
" [2 , 1] = -1.0\n",
" [501 , 1] = -1.0\n",
" [1 , 2] = -1.0\n",
" [2 , 2] = 4.0\n",
" [3 , 2] = -1.0\n",
" [502 , 2] = -1.0\n",
" [2 , 3] = -1.0\n",
" [3 , 3] = 4.0\n",
" [4 , 3] = -1.0\n",
" [503 , 3] = -1.0\n",
" [3 , 4] = -1.0\n",
" ⋮\n",
" [249997, 249997] = 4.0\n",
" [249998, 249997] = -1.0\n",
" [249498, 249998] = -1.0\n",
" [249997, 249998] = -1.0\n",
" [249998, 249998] = 4.0\n",
" [249999, 249998] = -1.0\n",
" [249499, 249999] = -1.0\n",
" [249998, 249999] = -1.0\n",
" [249999, 249999] = 4.0\n",
" [250000, 249999] = -1.0\n",
" [249500, 250000] = -1.0\n",
" [249999, 250000] = -1.0\n",
" [250000, 250000] = 4.0"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SparseArrays\n",
"A = spdiagm(-M => a_S[:][(1+M):end ],\n",
" -1 => a_W[:][ 2:end ],\n",
" 0 => a_P[:],\n",
" 1 => a_E[:][ 1:end-1],\n",
" M => a_N[:][ 1:end-M])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`SparseMatrixCSC` is a (parametric) composite type which is employed to store efficiently only the nonzero elements of a sparse matrix:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(:m, :n, :colptr, :rowval, :nzval)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fieldnames(SparseMatrixCSC)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"250000"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A.m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sparse direct solvers are also defined in `LinearAlgebra` module, therefore we can solve the considered linear system by simply using the usual backslash `\\` operator. In this particular problem, for simplicity, we are not enforcing a hot spot at a `hot_spot_location` but we are defining a heat source at the same location, namely `heat_source_location`:"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.971391 seconds (54 allocations: 268.415 MiB, 0.36% gc time)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
"0.0 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"1.0 \n",
" \n",
"\n",
" \n",
"\n",
" \n",
"\n",
"0 \n",
" \n",
"\n",
"0.1 \n",
" \n",
"\n",
"0.2 \n",
" \n",
"\n",
"0.3 \n",
" \n",
"\n",
"0.4 \n",
" \n",
"\n",
"0.5 \n",
" \n",
"\n",
"0.6 \n",
" \n",
"\n",
"0.7 \n",
" \n",
"\n",
"0.8 \n",
" \n",
"\n",
"0.9 \n",
" \n",
"\n",
"1.0 \n",
" \n",
" \n"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"source = zeros(M, M)\n",
"heat_source_location = (125, 100:400)\n",
"source[heat_source_location...] .= 1\n",
"@time T = A\\source[:]\n",
"\n",
"x = y = range(0, 1, length = M)\n",
"h3 = heatmap(x, y, reshape(T, M, M) / maximum(T),\n",
" xlabel = \"\\$x\\$\", xlim = (0,1),\n",
" ylabel = \"\\$y\\$\", ylim = (0,1),\n",
" tickfont = (:times,10), aspect_ratio = 1, clim=(0,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that $\\mathbf{A}$ is symmetric positive definite (SPD). The backslash operator `\\` recognizes this typical structure in order to choose the most suitable solver. If we change the sign of $\\mathbf{A}$ the solver will take longer since the matrix is not SPD any more:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.829143 seconds (54 allocations: 268.415 MiB, 0.28% gc time)\n",
" 1.564378 seconds (71 allocations: 402.102 MiB, 0.06% gc time)\n",
" 3.185904 seconds (71 allocations: 452.507 MiB, 0.05% gc time)\n"
]
}
],
"source": [
"@time T = A\\source[:] ;\n",
"A1 = (-1)*A\n",
"A2 = -A\n",
"@time T = A1\\source[:] ;\n",
"@time T = A2\\source[:] ;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rzamolo@units.it"
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
![\"drawing\"](\"FDgrid.png\")
\n",
"\n",
"In the FD method the derivatives in the heat equation are approximated as follows, where the subscript $_{i,j}$ is used for the spatial indexing on the cartesian lattice and the superscript $^{(n)}$ is used for the time indexing, where the time is $t=n\\Delta t$ with time step $\\Delta t$:\n",
"\n",
"$\\qquad\\displaystyle\\frac{\\partial \\phi}{\\partial t}\\Big\\vert_{i,j}^{(n)}\\approx\\displaystyle\\frac{\\phi_{i,j}^{(n+1)}-\\phi_{i,j}^{(n)}}{\\Delta t}$\n",
"\n",
"$\\qquad\\displaystyle\\Big(\\frac{\\partial^2\\phi}{\\partial x^2}+\\frac{\\partial^2\\phi}{\\partial y^2}\\Big)\\Big\\vert_{i,j}^{(n)}\\approx\\displaystyle\\frac{\\phi_{i+1,j}^{(n)}-2\\phi_{i,j}^{(n)}+\\phi_{i-1,j}^{(n)}}{\\Delta \n",
"x^2}+\\frac{\\phi_{i,j+1}^{(n)}-2\\phi_{i,j}^{(n)}+\\phi_{i,j-1}^{(n)}}{\\Delta \n",
"y^2}$\n",
"\n",
"By solving the approximated heat equation for $\\phi_{i,j}^{(n+1)}$, i.e., the solution at each node $i,j$ at the new time step $(n+1)$, we get the following scheme to perform an explicit time integration of the heat equation starting from an initial solution $\\phi(x,y,t=0)$:\n",
"\n",
"$\\phi_{i,j}^{(n+1)}=\\phi_{i,j}^{(n)}+\\Delta t \\displaystyle\\frac{\\phi_{i+1,j}^{(n)}-2\\phi_{i,j}^{(n)}+\\phi_{i-1,j}^{(n)}}{\\Delta \n",
"x^2}+\\Delta t\\frac{\\phi_{i,j+1}^{(n)}-2\\phi_{i,j}^{(n)}+\\phi_{i,j-1}^{(n)}}{\\Delta \n",
"y^2}$\n",
"\n",
"Boundary conditions must be enforced at the boundary $\\partial\\Omega$: for simplicity we consider homogeneous Dirichlet BCs, i.e., $\\phi_{i,j}=0$ for $i,j\\in\\{0,N\\}$.\n",
"\n",
"\n",
"The following scheme can be implemented very easily in Julia. Let's initialize $\\phi$ with (positive) random data for example:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 100\n",
"φ = zeros(N, N)\n",
"i = j = 2:(N-1)\n",
"φ[i,j] .= rand(N-2, N-2)\n",
"\n",
"Δx = 1 / (N-1)\n",
"Δt = Δx^2 / 5\n",
"for n ∈ 0:500\n",
" φ[i,j] += Δt * ( φ[i.+1,j] + φ[i.-1,j] + φ[i,j.+1] + φ[i,j.-1] - 4 * φ[i,j] ) / Δx^2\n",
"end\n",
"\n",
"x = y = range(0, 1, length = N)\n",
"heatmap(x, y, φ, xlim = (0,1), ylim = (0, 1), aspect_ratio = 1, tickfont = (:times,10), clim=(0,1)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If $\\phi$ is interpreted as a temperature field, we are solving an unsteady heat conduction problem on a square domain whose edges are constantly kept \"cool\" at $\\phi=0$. We can enforce a hot spot inside $\\Omega$ (e.g., constant \"hot\" temperature $\\phi=1$) and create an animation to figure out which is the time-dependent temperature distribution during the heating:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"