{
"cells": [
{
"cell_type": "markdown",
"id": "8863166a-db90-444f-af0a-e35b850c5492",
"metadata": {},
"source": [
"https://discourse.julialang.org/t/ode-solvers-why-is-matlab-ode45-uncannily-stable/63052/15"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "7dcd59cf-162d-4f64-bcd1-aead30056d89",
"metadata": {},
"outputs": [],
"source": [
"using DifferentialEquations\n",
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "b8917ff2-e1d2-467d-8cea-af15bbbca3ed",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"wk5! (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function wk5!(dP, P, params, t)\n",
"\n",
" #=\n",
"\n",
" The 5-Element WK with serial L\n",
" set Ls = 0 for 4-Element WK parallel\n",
"\n",
" Formulation for DifferentialEquations.jl\n",
"\n",
" P: solution vector (pressures p1 and p2)\n",
" params: parameter tuple\n",
" (Rc, Rp, C, Lp, Ls, I, q)\n",
"\n",
"\n",
" I need to find a way to tranfer the function name as well\n",
" for the time being we have to have the function in \"I\"\n",
"\n",
" =#\n",
"\n",
" # Split parameter tuple:\n",
" Rc, Rp, C, Lp, Ls, I, q = params\n",
"\n",
" dP[1] = (\n",
" -Rc / Lp * P[1]\n",
" + (Rc / Lp - 1 / Rp / C) * P[2]\n",
" + Rc * (1 + Ls / Lp) * didt(I, t, q)\n",
" + I(t, q) / C\n",
" )\n",
"\n",
" dP[2] = -1 / Rp / C * P[2] + I(t, q) / C\n",
"\n",
" return\n",
"\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "77c6203c-c7d9-4ecb-9bbe-e85d64d482aa",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"I_generic (generic function with 1 method)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Generic Input Waveform\n",
"# max volume flow in ml/s\n",
"max_i = 425\n",
"\n",
"# min volume flow in m^3/s\n",
"min_i = 0.0\n",
"\n",
"T = 0.9\n",
"\n",
"# Syst. Time in s\n",
"systTime = 2 / 5 * T\n",
"\n",
"# Dicrotic notch time\n",
"dicrTime = 0.02\n",
"\n",
"q_generic = (max_i, min_i, T, systTime, dicrTime)\n",
"\n",
"function I_generic(t, q_generic)\n",
" max_i, min_i, T, systTime, dicrTime = q_generic\n",
" # implicit conditional using boolean multiplicator\n",
" # sine waveform\n",
" (\n",
" (max_i - min_i) * sin(pi / systTime * (t % T))\n",
" * (t % T < (systTime + dicrTime) )\n",
" + min_i\n",
" )\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d6a69b16-c8da-4052-8c00-334ff59d92a7",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"didt (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function didt(I, t, q)\n",
" dt = 1e-3;\n",
" didt = (I(t+dt, q) - I(t-dt, q)) / (2 * dt)\n",
" return didt\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c4087c52-6218-48ed-9451-45f3f5331fc7",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(range(0, 2; length=2000), t -> didt(I_generic, t, q_generic); label=\"didt\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "40da522b-a3bf-4620-882d-f40602665824",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 29.388959 seconds (56.95 M allocations: 3.348 GiB, 2.89% gc time, 2.23% compilation time)\n",
" 1.840871 seconds (4.81 M allocations: 397.101 MiB, 8.46% gc time, 70.51% compilation time)\n",
" 2.746681 seconds (4.95 M allocations: 287.073 MiB, 99.96% compilation time)\n",
" 0.234024 seconds (1.34 M allocations: 133.246 MiB, 15.62% compilation time)\n",
" 2.391469 seconds (4.93 M allocations: 287.553 MiB, 99.98% compilation time)\n",
" 0.373095 seconds (1.94 M allocations: 192.746 MiB, 11.08% compilation time)\n",
" 11.395456 seconds (26.18 M allocations: 1.625 GiB, 2.22% gc time, 2.03% compilation time)\n"
]
}
],
"source": [
"# Initial condition and time span\n",
"P0 = [0.0, 0.0]\n",
"tspan = (0.0, 30.0)\n",
"\n",
"# Set parameters for Windkessel Model\n",
"Rc = 0.033\n",
"Rp = 0.6\n",
"C = 1.25\n",
"# L for serial model!\n",
"Ls = 0.01\n",
"# L for parallel\n",
"Lp = 0.02\n",
"\n",
"I = I_generic\n",
"q = q_generic\n",
"\n",
"p5 = (Rc, Rp, C, Lp, Ls, I, q)\n",
"\n",
"problem = ODEProblem(wk5!, P0, tspan, p5)\n",
"\n",
"dtmax = 1e-4\n",
"\n",
"@time solutionTsit = solve(problem); GC.gc()\n",
"@time solutionTsitLowdt = solve(problem, Tsit5(), dtmax=dtmax); GC.gc()\n",
"@time solutionBS3 = solve(problem, BS3()); GC.gc()\n",
"@time solutionBS3Lowdt = solve(problem, BS3(), dtmax=dtmax); GC.gc()\n",
"@time solutionDP5 = solve(problem,DP5()); GC.gc()\n",
"@time solutionDP5Lowdt = solve(problem,DP5(), dtmax=dtmax); GC.gc()\n",
"@time solutionStiff = solve(problem, alg_hints=[:stiff]); GC.gc()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4548ec77-bbb8-4154-a8f2-04ed6a6de0a0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.000408 seconds (1.88 k allocations: 195.938 KiB)\n",
" 0.349839 seconds (2.70 M allocations: 275.003 MiB)\n",
" 0.000939 seconds (3.96 k allocations: 444.500 KiB)\n",
" 0.186747 seconds (1.20 M allocations: 123.943 MiB)\n",
" 0.000301 seconds (1.09 k allocations: 113.859 KiB)\n",
" 0.325129 seconds (1.80 M allocations: 183.442 MiB)\n",
" 0.016856 seconds (104.41 k allocations: 4.175 MiB)\n"
]
}
],
"source": [
"@time solutionTsit = solve(problem); GC.gc()\n",
"@time solutionTsitLowdt = solve(problem, Tsit5(), dtmax=dtmax); GC.gc()\n",
"@time solutionBS3 = solve(problem, BS3()); GC.gc()\n",
"@time solutionBS3Lowdt = solve(problem, BS3(), dtmax=dtmax); GC.gc()\n",
"@time solutionDP5 = solve(problem,DP5()); GC.gc()\n",
"@time solutionDP5Lowdt = solve(problem,DP5(), dtmax=dtmax); GC.gc()\n",
"@time solutionStiff = solve(problem, alg_hints=[:stiff]); GC.gc()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "05cd24ba-5d04-4df6-8aad-8ca68648b3dd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"length(solutionTsit.t) = 192\n",
"length(solutionTsitLowdt.t) = 300010\n",
"length(solutionBS3.t) = 976\n",
"length(solutionBS3Lowdt.t) = 300034\n",
"length(solutionDP5.t) = 173\n",
"length(solutionDP5Lowdt.t) = 300002\n",
"length(solutionStiff.t) = 3992\n"
]
}
],
"source": [
"@show length(solutionTsit.t)\n",
"@show length(solutionTsitLowdt.t)\n",
"@show length(solutionBS3.t)\n",
"@show length(solutionBS3Lowdt.t)\n",
"@show length(solutionDP5.t)\n",
"@show length(solutionDP5Lowdt.t)\n",
"@show length(solutionStiff.t);"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "3373d20e-ae6f-446e-a7a4-2af529dcf793",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a, b = 0, 2\n",
"\n",
"plot()\n",
"#plot!(t -> solutionTsit(t; idxs=1), a, b; label=\"Tsit\")\n",
"plot!(t -> solutionTsitLowdt(t; idxs=1), a, b; label=\"TsitLowdt\")\n",
"#plot!(t -> solutionBS3(t; idxs=1), a, b; label=\"BS3\", ls=:dash)\n",
"plot!(t -> solutionBS3Lowdt(t; idxs=1), a, b; label=\"BS3Lowdt\", ls=:dash)\n",
"#plot!(t -> solutionDP5(t; idxs=1), a, b; label=\"DP5\", ls=:dashdot)\n",
"plot!(t -> solutionDP5Lowdt(t; idxs=1), a, b; label=\"DP5Lowdt\", ls=:dashdot)\n",
"plot!(t -> solutionStiff(t; idxs=1), a, b; label=\"Stiff\", ls=:dot, lw=1.5)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "3506abe4-cd40-4086-a148-8838ea114e19",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a, b = 0, 2\n",
"\n",
"plot()\n",
"plot!(t -> solutionTsit(t; idxs=1), a, b; label=\"Tsit\")\n",
"#plot!(t -> solutionTsitLowdt(t; idxs=1), a, b; label=\"TsitLowdt\")\n",
"plot!(t -> solutionBS3(t; idxs=1), a, b; label=\"BS3\", ls=:dash)\n",
"#plot!(t -> solutionBS3Lowdt(t; idxs=1), a, b; label=\"BS3Lowdt\", ls=:dash)\n",
"plot!(t -> solutionDP5(t; idxs=1), a, b; label=\"DP5\", ls=:dashdot)\n",
"#plot!(t -> solutionDP5Lowdt(t; idxs=1), a, b; label=\"DP5Lowdt\", ls=:dashdot)\n",
"plot!(t -> solutionStiff(t; idxs=1), a, b; label=\"Stiff\", ls=:dot, lw=1.5)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "fbeed588-52b2-4d23-90dd-cfd15ba19ce9",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a, b = 28, 30\n",
"\n",
"plot()\n",
"plot!(t -> solutionTsit(t; idxs=1), a, b; label=\"Tsit\")\n",
"#plot!(t -> solutionTsitLowdt(t; idxs=1), a, b; label=\"TsitLowdt\")\n",
"plot!(t -> solutionBS3(t; idxs=1), a, b; label=\"BS3\", ls=:dash)\n",
"#plot!(t -> solutionBS3Lowdt(t; idxs=1), a, b; label=\"BS3Lowdt\", ls=:dash)\n",
"plot!(t -> solutionDP5(t; idxs=1), a, b; label=\"DP5\", ls=:dashdot)\n",
"#plot!(t -> solutionDP5Lowdt(t; idxs=1), a, b; label=\"DP5Lowdt\", ls=:dashdot)\n",
"plot!(t -> solutionStiff(t; idxs=1), a, b; label=\"Stiff\", ls=:dot, lw=1.5)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a7f02bcd-d6b2-47aa-86c8-82bf9f45bbae",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"jupytext": {
"formats": "ipynb,auto:hydrogen"
},
"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"
}
},
"nbformat": 4,
"nbformat_minor": 5
}