Skip to content

Instantly share code, notes, and snippets.

@shivin9

shivin9/dextorious_Neumann.ipynb

Last active June 28, 2017 06:27
Show Gist options
  • Save shivin9/cd94aac140f9cf13183d88db334cc163 to your computer and use it in GitHub Desktop.
Save shivin9/cd94aac140f9cf13183d88db334cc163 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
dextorious_Neumann.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
dextorious_Neumann_aorder=3.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"using BenchmarkTools, DifferentialEquations, Plots\n",
"pyplot();"
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"heateqsystem (generic function with 1 method)"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function heateqsystem(t, u, du)\n",
" du[1] = -4.722*u[1] + 6.0*u[2] - 1.5*u[3] + 0.222*u[4] - 3.66*1/length(u)\n",
" du[2] = +0.9167*u[1] - 1.667*u[2] + 0.5*u[3] + 0.33 *u[4] - 0.0833*u[5]\n",
" du[3] = -0.0833*u[1] - 1.333*u[2] - 2.5*u[3] + 1.33 *u[4] - 0.0833*u[5]\n",
" du[4] = -0.0833*u[1] + 0.333*u[2] + 0.5*u[3] - 1.667 *u[4] + 0.91667*u[5]\n",
"\n",
" @fastmath @inbounds for i in 5 : length(u)-4\n",
" du[i] = -0.833*u[i-2] + 1.33*u[i+1] - 2.5*u[i] + 1.33*u[i-1] - 0.833*u[i-2]\n",
" end\n",
" \n",
" du[end-3] = -0.0833*u[end-4] + 1.333*u[end-3] - 2.5*u[end-2] + 1.33 *u[end-1] - 0.0833*u[end]\n",
" du[end-2] = -0.0833*u[end-4] + 1.333*u[end-3] - 2.5*u[end-2] + 1.33 *u[end-1] - 0.0833*u[end]\n",
" du[end-1] = -0.0833*u[end-4] + 1.333*u[end-3] - 2.5*u[end-2] + 1.33 *u[end-1] - 0.0833*u[end]\n",
" du[end] = -4.722*u[end-3] + 6.0*u[end-2] - 1.5*u[end-1] + 0.222*u[end] - 3.66*1/length(u)\n",
"\n",
" scale!(du, length(u)^2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"test (generic function with 1 method)"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function test(N)\n",
" h = 1 / (N-1)\n",
" x = collect(0 : h : 1)\n",
" u0 = -(x - 0.5).^2 + 1/12\n",
" tspan = (0., 1.)\n",
" prob1 = ODEProblem(heateqsystem, u0, tspan)\n",
" sol1 = solve(prob1, dense=false, saveat=[0.03, 1.0], dt=0.01,adaptive=false)\n",
" return x, sol1, prob1\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mInstability detected. Aborting\u001b[39m\n"
]
}
],
"source": [
"x, sol1, prob1 = test(100);"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"ename": "LoadError",
"evalue": "\u001b[91mSolution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.\u001b[39m",
"output_type": "error",
"traceback": [
"\u001b[91mSolution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.\u001b[39m",
"",
"Stacktrace:",
" [1] \u001b[1mode_interpolation\u001b[22m\u001b[22m at \u001b[1m/home/shivin/.julia/v0.6/OrdinaryDiffEq/src/dense/generic_dense.jl:225\u001b[22m\u001b[22m [inlined]",
" [2] \u001b[1m(::OrdinaryDiffEq.InterpolationData{#heateqsystem,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}})\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Float64, ::Void, ::Type{T} where T\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/shivin/.julia/v0.6/OrdinaryDiffEq/src/interp_func.jl:57\u001b[22m\u001b[22m",
" [3] \u001b[1m(::DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#heateqsystem,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#heateqsystem,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}})\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Float64\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/shivin/.julia/v0.6/DiffEqBase/src/solutions/ode_solutions.jl:14\u001b[22m\u001b[22m"
]
}
],
"source": [
"plot(x, sol1(0.5))"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100-element Array{Any,1}:\n",
" 1.0101\n",
" -14.5513\n",
" -30.1127\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" ⋮ \n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072\n",
" -45.2072"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"first_order_coeffs = [-2.083333, 4.0, -3.0, 1.333333, -0.25]\n",
"derivs_start = []\n",
"derivs_end = []\n",
"\n",
"for i in 0:1/99:1\n",
" push!(derivs_start, sum(first_order_coeffs.*sol1(i)[1:5])*100)\n",
" push!(derivs_end, sum(first_order_coeffs.*sol1(i)[end-4:end])*100)\n",
"end\n",
"derivs_start\n",
"# plot(x,derivs_start)\n",
"# plot(x, derivs_end)\n",
"# derivs_end"
]
}
],
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment