{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lineáris algebra"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:54.268000+02:00",
"start_time": "2019-04-09T08:24:54.226Z"
},
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Int64,2}:\n",
" 3 2 3\n",
" 1 4 1\n",
" 4 4 2"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Random mátrix\n",
"A = rand(1:4,3,3)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:54.500000+02:00",
"start_time": "2019-04-09T08:24:54.491Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Egyesekkel teli vektor\n",
"x = fill(1.0, (3)) # vagy `ones(3)`"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:55.083000+02:00",
"start_time": "2019-04-09T08:24:55.028Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 8.0\n",
" 6.0\n",
" 10.0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = A*x # Szorzás"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:56.498000+02:00",
"start_time": "2019-04-09T08:24:55.534Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Adjoint{Int64,Array{Int64,2}}:\n",
" 3 1 4\n",
" 2 4 4\n",
" 3 1 2"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Konjugált transzponált\n",
"A'"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:57.607000+02:00",
"start_time": "2019-04-09T08:24:56.639Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Transpose{Int64,Array{Int64,2}}:\n",
" 3 1 4\n",
" 2 4 4\n",
" 3 1 2"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Transzponált\n",
"transpose(A)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:58.836000+02:00",
"start_time": "2019-04-09T08:24:57.594Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 1.0000000000000002\n",
" 1.0 \n",
" 0.9999999999999998"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Lineáris egyenletrendszerek megoldása\n",
"A\\b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mátrix dekompozíciók"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:58.837000+02:00",
"start_time": "2019-04-09T08:24:58.402Z"
}
},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.268823 0.186721 0.321859\n",
" 0.450096 0.186799 0.961101\n",
" 0.134243 0.964543 0.176886"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Random mátrix\n",
"A = rand(3,3)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.2890541723184847"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigmax(A)"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}\n",
"eigenvalues:\n",
"3-element Array{Float64,1}:\n",
" 1.2890541723184847 \n",
" 0.10173046819308784\n",
" -0.7582765663881217 \n",
"eigenvectors:\n",
"3×3 Array{Float64,2}:\n",
" -0.32913 -0.904794 -0.101461\n",
" -0.69335 0.0935566 -0.685181\n",
" -0.641045 0.415446 0.721272"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigen(A)"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:58.837000+02:00",
"start_time": "2019-04-09T08:24:58.402Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}\n",
"Q factor:\n",
"3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:\n",
" -0.496736 -0.0634131 -0.865582 \n",
" -0.831697 -0.250269 0.495626 \n",
" -0.248058 0.966097 0.0715773\n",
"R factor:\n",
"3×3 Array{Float64,2}:\n",
" -0.541178 -0.487374 -1.0031 \n",
" 0.0 0.873252 -0.0900548\n",
" 0.0 0.0 0.210412 "
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"qr(A)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"LU{Float64,Array{Float64,2}}\n",
"L factor:\n",
"3×3 Array{Float64,2}:\n",
" 1.0 0.0 0.0\n",
" 0.298255 1.0 0.0\n",
" 0.597256 0.0826928 1.0\n",
"U factor:\n",
"3×3 Array{Float64,2}:\n",
" 0.450096 0.186799 0.961101\n",
" 0.0 0.908829 -0.109767\n",
" 0.0 0.0 -0.243087"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lu(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Speciális Mátrix struktúrák"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:58.837000+02:00",
"start_time": "2019-04-09T08:24:58.402Z"
}
},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:24:58.940000+02:00",
"start_time": "2019-04-09T08:24:58.866Z"
}
},
"outputs": [],
"source": [
"n = 1000\n",
"A = randn(n,n);"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:25:00.362000+02:00",
"start_time": "2019-04-09T08:25:00.177Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Szimmetria ellenőrzése\n",
"Asym = A + A'\n",
"issymmetric(Asym)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:25:00.931000+02:00",
"start_time": "2019-04-09T08:25:00.890Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"false"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Numerikus zajjal terhelt mátrix\n",
"Asym_noisy = copy(Asym)\n",
"Asym_noisy[1,2] += 5eps()\n",
"issymmetric(Asym_noisy)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:25:01.286000+02:00",
"start_time": "2019-04-09T08:25:01.269Z"
}
},
"outputs": [],
"source": [
"# Szimmetrikus mátrixok direkt definiálása\n",
"Asym_explicit = Symmetric(Asym_noisy);"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:25:03.276000+02:00",
"start_time": "2019-04-09T08:25:01.585Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.090288 seconds (18.99 k allocations: 8.979 MiB)\n",
" 0.862438 seconds (18 allocations: 7.921 MiB)\n",
" 0.197209 seconds (7.60 k allocations: 8.378 MiB)\n"
]
}
],
"source": [
"@time eigvals(Asym);\n",
"@time eigvals(Asym_noisy);\n",
"@time eigvals(Asym_explicit);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Nagy elemszámú problémák"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:25:08.334000+02:00",
"start_time": "2019-04-09T08:25:07.649Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.873977 seconds (519.48 k allocations: 208.403 MiB, 8.45% gc time)\n"
]
},
{
"data": {
"text/plain": [
"6.334477063028135"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 1_000_000;\n",
"A = SymTridiagonal(randn(n), randn(n-1));\n",
"@time eigmax(A)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"ExecuteTime": {
"end_time": "2019-04-09T10:25:38.049000+02:00",
"start_time": "2019-04-09T08:25:38.044Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"1000000×1000000 SymTridiagonal{Float64,Array{Float64,1}}:\n",
" 0.693507 -0.511006 ⋅ ⋅ … ⋅ ⋅ ⋅ \n",
" -0.511006 -0.0556827 0.342219 ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ 0.342219 0.78431 -1.09116 ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ -1.09116 -0.62924 ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ 1.52308 ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋮ ⋱ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ 0.458042 ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ 0.559203 1.33234 ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ 1.33234 -0.81644 0.63982 \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ 0.63982 0.0101003"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Differenciálegyenletek numerikus megoldása\n",
"\n",
"[Differentialequations.jl](http://docs.juliadiffeq.org/latest/) - Beállítások, és részletek a megoldókról általában"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"using DifferentialEquations\n",
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\n",
"# Probléma definiálása (itt most Közönséges DE kezdeti érték problémája)\n",
"f(u,p,t) = 1.01*u\n",
"u0=1/2\n",
"tspan = (0.0,1.0)\n",
"prob = ODEProblem(f,u0,tspan)\n",
"\n",
"# Megoldása\n",
"sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)\n",
"\n",
"\n",
"# Ábrázolása\n",
"plot(sol,linewidth=5,title=\"Solution to the linear ODE with a thick line\",\n",
" xaxis=\"Time (t)\",yaxis=\"u(t) (in μm)\",label=\"My Thick Line!\") # legend=false\n",
"plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label=\"True Solution!\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lorentz egyenletek\n",
"\n",
"\\begin{align}\n",
"\\frac{dx}{dt} &= σ(y-x) \\\\\n",
"\\frac{dy}{dt} &= x(ρ-z) - y \\\\\n",
"\\frac{dz}{dt} &= xy - βz \\\\\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"ename": "UndefVarError",
"evalue": "UndefVarError: parameterized_lorenz not defined",
"output_type": "error",
"traceback": [
"UndefVarError: parameterized_lorenz not defined",
"",
"Stacktrace:",
" [1] top-level scope at In[23]:4"
]
}
],
"source": [
"u0 = [1.0,0.0,0.0]\n",
"tspan = (0.0,1.0)\n",
"p = [10.0,28.0,8/3]\n",
"prob = ODEProblem(parameterized_lorenz,u0,tspan,p)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lorenz (generic function with 1 method)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function lorenz(du,u,p,t)\n",
" du[1] = p[1]*(u[2]-u[1])\n",
" du[2] = u[1]*(p[2]-u[3]) - u[2]\n",
" du[3] = u[1]*u[2] - (p[3])*u[3]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 6.783799 seconds (17.47 M allocations: 899.842 MiB, 9.03% gc time)\n"
]
}
],
"source": [
"u0 = [1.0;0.0;0.0]\n",
"tspan = (0.0,100.0)\n",
"p=(15.,28.,8. /3.)\n",
"prob = ODEProblem(lorenz,u0,tspan,p)\n",
"@time sol = solve(prob);"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol,vars=(1,2,3))"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol,vars=(0,2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Nemlineáris egyenletek megoldása"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Általánosabb nemlineáris megoldó\n",
"\n",
"[(NLsolve.jl)](https://github.com/JuliaNLSolvers/NLsolve.jl)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"using NLsolve"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Results of Nonlinear Solver Algorithm\n",
" * Algorithm: Trust-region with dogleg and autoscaling\n",
" * Starting Point: [0.1, 1.2]\n",
" * Zero: [-3.7818e-16, 1.0]\n",
" * Inf-norm of residuals: 0.000000\n",
" * Iterations: 4\n",
" * Convergence: true\n",
" * |x - x'| < 0.0e+00: false\n",
" * |f(x)| < 1.0e-08: true\n",
" * Function Calls (f): 5\n",
" * Jacobian Calls (df/dx): 5"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Probléma\n",
"function f!(F, x)\n",
" F[1] = (x[1]+3)*(x[2]^3-7)+18\n",
" F[2] = sin(x[2]*exp(x[1])-1)\n",
"end\n",
"\n",
"# Probléma Jacobi mátrixa\n",
"function j!(J, x)\n",
" J[1, 1] = x[2]^3-7\n",
" J[1, 2] = 3*x[2]^2*(x[1]+3)\n",
" u = exp(x[1])*cos(x[2]*exp(x[1])-1)\n",
" J[2, 1] = x[2]*u\n",
" J[2, 2] = u\n",
"end\n",
"# Kezdeti megoldásvektor\n",
"init=[ 0.1; 1.2];\n",
"solu=nlsolve(f!, j!, init)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Intervallum aritmetika segítségével\n",
"\n",
"[(IntervalRootFinding.jl)](https://github.com/JuliaIntervals/IntervalRootFinding.jl)"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"using IntervalArithmetic, IntervalRootFinding\n",
"using StaticArrays"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Root{Interval{Float64}},1}:\n",
" Root([1.73205, 1.73206], :unique) \n",
" Root([1.41421, 1.41422], :unknown) \n",
" Root([-1.41422, -1.41421], :unknown)\n",
" Root([-1.73206, -1.73205], :unique) "
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g(x) = (x^2-2)^2 * (x^2 - 3)\n",
"rts=roots(g, -10..10)"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.027744 seconds (4.80 M allocations: 224.263 MiB, 7.49% gc time)\n"
]
},
{
"data": {
"text/plain": [
"4-element Array{Root{IntervalBox{2,Float64}},1}:\n",
" Root([1.94053, 1.94054] × [1.49727, 1.49728], :unique) \n",
" Root([1.60901, 1.60902] × [1.45725, 1.45726], :unique) \n",
" Root([1.10116, 1.10117] × [1.377, 1.37701], :unique) \n",
" Root([-2.00126e-16, 4.06797e-16] × [0.999999, 1.00001], :unique)"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function g(x)\n",
" (x1, x2) = x\n",
" SVector((x1+3)*(x2^3-7)+18,\n",
" sin(x2*exp(x1)-1))\n",
"end\n",
"interval=-2..2\n",
"@time rts = roots(g, interval × interval) # × = \\times"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Felezőmódszer (Multi Dimensional Bisection Method)\n",
"\n",
"https://github.com/bachrathyd/MDBM.jl\n",
"\n",
"Megjegyzés: Ez a módszer inkább ponthalmazok megtalálására alkalmas, pl stabilitási határok"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"using MDBM"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Array{Float64,1},1}:\n",
" [1.5842e-5, -1.56338e-5, -1.56143e-5, 1.53927e-5, 1.1012, 1.10117, 1.10113, 1.10118, 1.60902, 1.60905, 1.94062, 1.94049, 1.94053, 1.9405]\n",
" [1.00002, 0.999977, 0.999976, 1.00002, 1.37705, 1.37698, 1.37701, 1.377, 1.45723, 1.45723, 1.49728, 1.49729, 1.49726, 1.49731] "
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function foo(x1,x2)\n",
" SVector((x1+3)*(x2^3-7)+18,\n",
" sin(x2*exp(x1)-1))\n",
"end\n",
"\n",
"ax1=Axis(-2.0:0.1:2.0) # initial grid in x1 direction\n",
"ax2=Axis(-2.0:0.1:2.0) # initial grid in x2 direction\n",
"\n",
"mymdbm=MDBM_Problem(foo,[ax1,ax2])\n",
"iteration=3 #number of refinements (resolution doubling)\n",
"MDBM.solve!(mymdbm,iteration)\n",
"\n",
"getinterpolatedsolution(mymdbm)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter(getinterpolatedsolution(mymdbm)...,size=(300,200),label=\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Késleltetett egyenletek stabilitásvizsgálata\n",
"\n",
"https://github.com/HTSykora/SemiDiscretizationMethod.jl"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"using SemiDiscretizationMethod\n",
"using MDBM\n",
"using Plots\n",
"using LaTeXStrings"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Késleltetett Mathieu egyenlet\n",
"\n",
"$$\\ddot{x}(t) + a_1 \\,\\dot{x}(t)+(\\delta + \\varepsilon \\cos(t))x(t)= b_0 \\,x(t-2\\pi) + \\sin(2t)$$\n",
"Elsőrendű formára átírva:
\n",
"$$\\dot{\\mathbf{x}}(t) = \\mathbf{A}(t) \\mathbf{x}(t) + \\sum_{j=1}^g \\mathbf{B}(t) \\mathbf{x}(t-\\tau_j(t))+\\mathbf{c}(t)$$\n",
"
\n",
"$$ \\mathbf{x}(t) = \\begin{bmatrix}x(t) \\\\ \\dot{x}(t)\\end{bmatrix}, \\quad\n",
"\\mathbf{A}(t) = \\begin{bmatrix} 0 & 1 \\\\ -\\delta - \\varepsilon \\cos(t) & -a_1 \\end{bmatrix},\n",
"\\quad \\mathbf{B}_1(t) = \\begin{bmatrix}0 & 0 \\\\ b_0 & 0\\end{bmatrix},\n",
"\\quad \\tau_1(t) \\equiv 2\\pi, \n",
"\\quad \\text{and} \\quad \\mathbf{c}(t) = \\begin{bmatrix} 0 \\\\ \\sin(2t) \\end{bmatrix}.$$"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"function createMathieuProblem(δ,ε,b0,a1;T=2π)\n",
" AMx = ProportionalMX(t->@SMatrix [0. 1.; -δ-ε*cos(2π/T*t) -a1]);\n",
" τ1=2π # if function is needed, the use τ1 = t->foo(t)\n",
" BMx1 = DelayMX(τ1,t->@SMatrix [0. 0.; b0 0.]);\n",
" cVec = Additive(t->@SVector [0.,sin(4π/T*t)])\n",
" LDDEProblem(AMx,[BMx1],cVec)\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"τmax=2π # the largest τ of the system\n",
"P=2π #Principle period of the system (sin(t)=sin(t+P)) \n",
"mathieu_lddep=createMathieuProblem(3.,2.,-0.15,0.1,T=P); # LDDE problem for Hayes equation\n",
"method=SemiDiscretization(1,0.01) # 3rd order semi discretization with Δt=0.1\n",
"# if P = τmax, then n_steps is automatically calculated\n",
"mapping=calculateMapping(mathieu_lddep,method,τmax,\n",
" n_steps=Int((P+100eps(P))÷method.Δt),calculate_additive=true); #The discrete mapping of the system\n",
"\n",
"a1=0.1;\n",
"ε=1;\n",
"τmax=2π;\n",
"T=1π;\n",
"method=SemiDiscretization(2,T/40);\n",
"\n",
"foo(δ,b0) = log(spectralRadiusOfMapping(calculateMapping(createMathieuProblem(δ,ε,b0,a1,T=T),method,τmax,\n",
" n_steps=Int((T+100eps(T))÷method.Δt)))); # No additive term calculated\n",
"\n",
"axis=[Axis(-1:0.2:5.,:δ),\n",
" Axis(-2:0.2:1.5,:b0)];"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 13.424935 seconds (43.76 M allocations: 11.089 GiB, 21.67% gc time)\n"
]
}
],
"source": [
"iteration=3;\n",
"@time stab_border_points=getinterpolatedsolution(MDBM.solve!(MDBM_Problem(foo,axis),iteration));"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"