Instantly share code, notes, and snippets.

# shivin9/robin_problems.ipynb

Last active July 14, 2017 11:34
Show Gist options
• Save shivin9/124ed1e5ea96792fc8666e0caf32715c to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 { "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Plots.PyPlotBackend()" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using DifferentialEquations, Plots, PDEOperators\n", "pyplot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Demonstrating the correctness of Neumann BC as it is used in Robin BC\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "N = 100;\n", "h_inv = 1/(N-1);\n", "x = collect(0 : h_inv : 1);\n", "u0 = -(x - 0.5).^2 + 1/12;\n", "plot(x,u0)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "left_deriv = sum(first_order_coeffs .* u0[1:4]) = 1.0000000000000107\n", "right_deriv = -(sum(reverse(first_order_coeffs) .* u0[end - 3:end])) = -0.9999999999999893\n" ] } ], "source": [ "first_order_coeffs = [-11/6, 3.0, -3/2, 1/3]*(N-1)\n", "@show left_deriv = sum(first_order_coeffs.*u0[1:4]);\n", "@show right_deriv = -sum(reverse(first_order_coeffs).*u0[end-3:end]); # should theoretically be -1\n", "A = LinearOperator{Float64}(2,2,1/(N-1),N,:Neumann,:Neumann;bndry_fn=(left_deriv, right_deriv));" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "prob1 = ODEProblem(A, u0, (0.,10.));\n", "sol1 = solve(prob1, CVODE_BDF(), dense=false, tstops=0:0.01:10);" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "derivs_start = []\n", "derivs_end = []\n", "NS = length(first_order_coeffs)\n", "t = 0:0.1:10\n", "for i in t\n", " push!(derivs_start, sum(first_order_coeffs.*sol1(i)[1 : NS]))\n", " push!(derivs_end, -sum(reverse(first_order_coeffs).*sol1(i)[end-NS+1 : end]))\n", "end\n", "plot(derivs_start)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(derivs_end)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Robin BC" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$u(1,t) + u'(1,t) = 0$$\n", "initiliized as $$u(0,0) + u'(0,0)$$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "A = LinearOperator{Float64}(2,2,1/(N-1),N,:Robin,:Dirichlet;bndry_fn=((1,1,-left_deriv+u0[1]),u0[end]));\n", "prob1 = ODEProblem(A, u0, (0.,10.));\n", "# sol1 = solve(prob1, CVODE_BDF(), dense=false, tstops=0:0.01:10); # giving instability error\n", "# sol_adams = solve(prob1, CVODE_Adams(), dense=false, tstops=0:0.01:10); # giving instability error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Strangely, it works after we change the domain to [$-\\pi$, $\\pi$]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x = collect(-pi : 2pi/1023 : pi);\n", "u0 = -(x - 0.5).^2 + 1/12;\n", "first_order_coeffs = [-11/6, 3.0, -3/2, 1/3] * (1023/2π)\n", "left_deriv = sum(first_order_coeffs.*u0[1:4])\n", "A = LinearOperator{Float64}(2,2,2π/1023,1024,:Robin,:Dirichlet;bndry_fn=((1,1,-left_deriv+u0[1]),u0[end]));" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "prob1 = ODEProblem(A, u0, (0.,10.));\n", "sol1 = solve(prob1, dense=false, tstops=0:0.01:10);\n", "# Strangely it doesn't work with CVODE_BDF() but does with CVODE_Adams()\n", "sol_CVODE = solve(prob1, CVODE_Adams(), dense=false, tstops=0:0.01:10);" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(x, [sol1(i) for i in 0:1:10])" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(x, [sol_CVODE(i) for i in 0:1:10])" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(x, [sol1(i) for i in 0:1:10].-[sol_CVODE(i) for i in 0:1:10])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "robin_normal = []\n", "robin_adams = []\n", "val_start = -left_deriv + u0[1]\n", "for i in 0:0.1:10\n", " push!(robin_normal, (-sum(first_order_coeffs.*sol1(i)[1:4])+sol1(i)[1]-val_start))\n", " push!(robin_adams, (-sum(first_order_coeffs.*sol_CVODE(i)[1:4])+sol1(i)[1]-val_start))\n", "end\n", "plot(robin_normal);\n", "plot!(robin_adams) # overlapping" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.0", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.0" } }, "nbformat": 4, "nbformat_minor": 2 }
to join this conversation on GitHub. Already have an account? Sign in to comment