linearization-of-single-variable-non-linear-function.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Linearization of Non-Linear Equations"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "using Plots, LinearAlgebra\n",
        "plotly()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Problem Statement\n",
        "\n",
        "Consider the following non-linear equation. \n",
        "\n",
        "$\\displaystyle{u\\left( x \\right) = x - x^2 + \\sqrt\\gamma\\left( \\frac{1}{4\\pi}\\sin\\left( \\frac{2\\pi x}{\\gamma} \\right)-\n",
        "\\frac{x}{2\\pi}\\sin\\left( \\frac{2\\pi x}{\\gamma} \\right) -\\frac{\\gamma}{4\\pi^2}\\cos\\left( \\frac{2\\pi x}{\\gamma} \\right) + \\frac{\\gamma}{4\\pi^2} \\right) }\n",
        "$\n",
        "\n",
        "where $\\gamma$ is a parameter that can take values from $10^{-1} \\text{ to } 10^{-2}$ and $x \\in [0,1]$."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "gm = 0.01:0.01:0.1;\n",
        "x_ini = 0.0:0.01:1.;\n",
        "γ = gm[1];\n",
        "x0 = x_ini[1]\n",
        "u(x, γ=γ) = x-x^2 +√γ*(1/4π * sin(2π*x/γ) - x/2π * sin(2π*x/γ) - γ/4π^2 * cos(2π*x/γ) + γ/4π^2);\n",
        "du(x) = (1-2x)*(1+ 1/(2*√γ) * cos(2π*x/γ));\n",
        "uₗ(x, x0=x0) = u(x0) + du(x0)*(x-x0);\n",
        "x = 0:1e-5:1;\n",
        "ϵ = 1e-11;\n",
        "x_star = 0.5;\n",
        "x_next(xi, γ=γ) = xi + (u(x_star, γ) - u(xi))/du(xi);"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "plot(u, x, legend=false)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "u_values = Array{Vector{Float64}}(undef,length(gm), length(x_ini))\n",
        "residue_values = Array{Vector{Float64}}(undef,length(gm), length(x_ini))\n",
        "x_values = Array{Vector{Float64}}(undef,length(gm), length(x_ini))\n",
        "ustar_values = Array{Float64}(undef,length(gm), length(x_ini))\n",
        "itr_values = Array{Int64}(undef, length(gm), length(x_ini))\n",
        "for i in 1:length(gm)\n",
        "    for j in 1:length(x_ini)\n",
        "        γ = gm[i];\n",
        "        x0 = x_ini[j];\n",
        "        xs = [x0]\n",
        "        xi = last(xs)\n",
        "        u_star = u(x_star, γ)\n",
        "        u_L = [uₗ(x_star, xi),]\n",
        "        residue = [abs(u_star - last(u_L))]\n",
        "        its = 0\n",
        "        while last(residue) > ϵ && its < 10^4\n",
        "            xi = x_next(xi, γ)\n",
        "            push!(xs, xi)\n",
        "            push!(u_L, uₗ(x_star, xi))\n",
        "            push!(residue, abs(u_star - last(u_L)))\n",
        "            its += 1\n",
        "        end\n",
        "        u_values[i, j] = u_L\n",
        "        residue_values[i,j] = residue\n",
        "        x_values[i,j] = xs\n",
        "        ustar_values[i,j] = u_star\n",
        "        itr_values[i,j] = length(residue)\n",
        "    end\n",
        "end"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "p1=surface(x_ini, gm, 1 ./itr_values, zlim=(0, 0.15), c=:seaborn_bright,# cgrad(:gnuplot, rev=true),\n",
        "        xlabel=\"x₀ value\", ylabel=\"γ value\", zlabel=\"1/Iterations\", cbar=false)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "contourf(x_ini, gm, itr_values, zlim=(0, 0.15), c=:seaborn_bright,# cgrad(:gnuplot, rev=true),\n",
        "        xlabel=\"x₀ value\", ylabel=\"γ value\")"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "i = rand(1:length(gm))\n",
        "j = rand(1:length(x_ini))\n",
        "itr_values[i,j]"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "u_star = ustar_values[i,j]\n",
        "xs = x_values[i,j]\n",
        "residue = residue_values[i,j]\n",
        "γ = gm[i]\n",
        "u_vals = u_values[i,j]\n",
        "p1=plot(x->u(x, γ), x, ylim=(0,0.3), xlabel=\"x\", ylabel=\"u(x)\")\n",
        "plot!(p1, x, fill(u_star, size(x)), lc=:red)\n",
        "c=0\n",
        "for xᵢ ∈ xs\n",
        "    plot!(p1, x->uₗ(x, xᵢ), x, ylim=(0,0.3), xlim=(0,1), legend=false, lc=:black)\n",
        "    plot!(p1, fill(x_next(xᵢ,γ), size(0:.05:.25)), 0:.05:.25, lc=:orange)\n",
        "    if length(xs) > 50\n",
        "        c += 1\n",
        "        if c > 10\n",
        "        break\n",
        "        end\n",
        "    end\n",
        "end\n",
        "p2 = plot(residue, scale=:log10, xlabel=\"Iteration\", ylabel=\"Residue\", legend=false)\n",
        "plot(p1, p2, layout=(2,1))\n",
        "p1"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "p2"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Julia 1.6.1",
      "language": "julia",
      "name": "julia-1.6"
    },
    "language_info": {
      "file_extension": ".jl",
      "mimetype": "application/julia",
      "name": "julia",
      "version": "1.6.1"
    },
    "nteract": {
      "version": "0.28.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}