{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Intégration numérique avec Julia, une introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"Il nous faut tout d'abord importer les bons packages de Julia"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Precompiling DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa]\n",
"└ @ Base loading.jl:1342\n",
"\u001b[91m\u001b[1mERROR: \u001b[22m\u001b[39mLoadError: ArgumentError: Package DataStructures does not have Compat in its dependencies:\n",
"- If you have DataStructures checked out for development and have\n",
" added Compat as a dependency but haven't updated your primary\n",
" environment's manifest file, try `Pkg.resolve()`.\n",
"- Otherwise you may need to report an issue with DataStructures\n",
"Stacktrace:\n",
" [1] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90minto\u001b[39m::\u001b[0mModule, \u001b[90mmod\u001b[39m::\u001b[0mSymbol\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:906\u001b[0m\n",
" [2] \u001b[0m\u001b[1minclude\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mBase.jl:386\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [3] \u001b[0m\u001b[1minclude_package_for_output\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90minput\u001b[39m::\u001b[0mString, \u001b[90mdepot_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mdl_load_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mload_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mconcrete_deps\u001b[39m::\u001b[0mVector\u001b[90m{Pair{Base.PkgId, UInt64}}\u001b[39m, \u001b[90msource\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1235\u001b[0m\n",
" [4] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
" [5] \u001b[0m\u001b[1meval\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mboot.jl:360\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [6] \u001b[0m\u001b[1meval\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mx\u001b[39m::\u001b[0mExpr\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase.MainInclude\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mclient.jl:446\u001b[0m\n",
" [7] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
"in expression starting at /Users/gergaud/.julia/packages/DataStructures/nBjdy/src/DataStructures.jl:1\n",
"\u001b[91m\u001b[1mERROR: \u001b[22m\u001b[39mLoadError: Failed to precompile DataStructures [864edb3b-99cc-5e75-8d2d-829cb0a9cfe8] to /Users/gergaud/.julia/compiled/v1.6/DataStructures/jl_FUK9uh.\n",
"Stacktrace:\n",
" [1] \u001b[0m\u001b[1merror\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90ms\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4merror.jl:33\u001b[0m\n",
" [2] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString, \u001b[90minternal_stderr\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90minternal_stdout\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90mignore_loaded_modules\u001b[39m::\u001b[0mBool\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1385\u001b[0m\n",
" [3] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1329\u001b[0m\n",
" [4] \u001b[0m\u001b[1m_require\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1043\u001b[0m\n",
" [5] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90muuidkey\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:936\u001b[0m\n",
" [6] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90minto\u001b[39m::\u001b[0mModule, \u001b[90mmod\u001b[39m::\u001b[0mSymbol\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:923\u001b[0m\n",
" [7] \u001b[0m\u001b[1minclude\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mBase.jl:386\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [8] \u001b[0m\u001b[1minclude_package_for_output\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90minput\u001b[39m::\u001b[0mString, \u001b[90mdepot_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mdl_load_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mload_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mconcrete_deps\u001b[39m::\u001b[0mVector\u001b[90m{Pair{Base.PkgId, UInt64}}\u001b[39m, \u001b[90msource\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1235\u001b[0m\n",
" [9] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
" [10] \u001b[0m\u001b[1meval\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mboot.jl:360\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [11] \u001b[0m\u001b[1meval\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mx\u001b[39m::\u001b[0mExpr\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase.MainInclude\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mclient.jl:446\u001b[0m\n",
" [12] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
"in expression starting at /Users/gergaud/.julia/packages/StatsBase/Q76Ni/src/StatsBase.jl:1\n",
"\u001b[91m\u001b[1mERROR: \u001b[22m\u001b[39mLoadError: Failed to precompile StatsBase [2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91] to /Users/gergaud/.julia/compiled/v1.6/StatsBase/jl_aTmMTx.\n",
"Stacktrace:\n",
" [1] \u001b[0m\u001b[1merror\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90ms\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4merror.jl:33\u001b[0m\n",
" [2] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString, \u001b[90minternal_stderr\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90minternal_stdout\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90mignore_loaded_modules\u001b[39m::\u001b[0mBool\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1385\u001b[0m\n",
" [3] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1329\u001b[0m\n",
" [4] \u001b[0m\u001b[1m_require\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1043\u001b[0m\n",
" [5] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90muuidkey\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:936\u001b[0m\n",
" [6] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90minto\u001b[39m::\u001b[0mModule, \u001b[90mmod\u001b[39m::\u001b[0mSymbol\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:923\u001b[0m\n",
" [7] \u001b[0m\u001b[1minclude\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mBase.jl:386\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [8] \u001b[0m\u001b[1minclude_package_for_output\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90minput\u001b[39m::\u001b[0mString, \u001b[90mdepot_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mdl_load_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mload_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mconcrete_deps\u001b[39m::\u001b[0mVector\u001b[90m{Pair{Base.PkgId, UInt64}}\u001b[39m, \u001b[90msource\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1235\u001b[0m\n",
" [9] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
" [10] \u001b[0m\u001b[1meval\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mboot.jl:360\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [11] \u001b[0m\u001b[1meval\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mx\u001b[39m::\u001b[0mExpr\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase.MainInclude\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mclient.jl:446\u001b[0m\n",
" [12] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
"in expression starting at /Users/gergaud/.julia/packages/Distributions/1cqVu/src/Distributions.jl:1\n",
"\u001b[91m\u001b[1mERROR: \u001b[22m\u001b[39mLoadError: Failed to precompile Distributions [31c24e10-a181-5473-b8eb-7969acd0382f] to /Users/gergaud/.julia/compiled/v1.6/Distributions/jl_tmg8Kh.\n",
"Stacktrace:\n",
" [1] \u001b[0m\u001b[1merror\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90ms\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4merror.jl:33\u001b[0m\n",
" [2] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString, \u001b[90minternal_stderr\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90minternal_stdout\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90mignore_loaded_modules\u001b[39m::\u001b[0mBool\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1385\u001b[0m\n",
" [3] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1329\u001b[0m\n",
" [4] \u001b[0m\u001b[1m_require\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1043\u001b[0m\n",
" [5] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90muuidkey\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:936\u001b[0m\n",
" [6] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90minto\u001b[39m::\u001b[0mModule, \u001b[90mmod\u001b[39m::\u001b[0mSymbol\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:923\u001b[0m\n",
" [7] \u001b[0m\u001b[1minclude\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mBase.jl:386\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [8] \u001b[0m\u001b[1minclude_package_for_output\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90minput\u001b[39m::\u001b[0mString, \u001b[90mdepot_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mdl_load_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mload_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mconcrete_deps\u001b[39m::\u001b[0mVector\u001b[90m{Pair{Base.PkgId, UInt64}}\u001b[39m, \u001b[90msource\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1235\u001b[0m\n",
" [9] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
" [10] \u001b[0m\u001b[1meval\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mboot.jl:360\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [11] \u001b[0m\u001b[1meval\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mx\u001b[39m::\u001b[0mExpr\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase.MainInclude\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mclient.jl:446\u001b[0m\n",
" [12] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
"in expression starting at /Users/gergaud/.julia/packages/DiffEqBase/b1nST/src/DiffEqBase.jl:1\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[91m\u001b[1mERROR: \u001b[22m\u001b[39mLoadError: Failed to precompile DiffEqBase [2b5f629d-d688-5b77-993f-72d75c75574e] to /Users/gergaud/.julia/compiled/v1.6/DiffEqBase/jl_L6uKwc.\n",
"Stacktrace:\n",
" [1] \u001b[0m\u001b[1merror\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90ms\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4merror.jl:33\u001b[0m\n",
" [2] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString, \u001b[90minternal_stderr\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90minternal_stdout\u001b[39m::\u001b[0mIOContext\u001b[90m{Base.PipeEndpoint}\u001b[39m, \u001b[90mignore_loaded_modules\u001b[39m::\u001b[0mBool\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1385\u001b[0m\n",
" [3] \u001b[0m\u001b[1mcompilecache\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90mpath\u001b[39m::\u001b[0mString\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1329\u001b[0m\n",
" [4] \u001b[0m\u001b[1m_require\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1043\u001b[0m\n",
" [5] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90muuidkey\u001b[39m::\u001b[0mBase.PkgId\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:936\u001b[0m\n",
" [6] \u001b[0m\u001b[1mrequire\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90minto\u001b[39m::\u001b[0mModule, \u001b[90mmod\u001b[39m::\u001b[0mSymbol\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:923\u001b[0m\n",
" [7] \u001b[0m\u001b[1minclude\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mBase.jl:386\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [8] \u001b[0m\u001b[1minclude_package_for_output\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mpkg\u001b[39m::\u001b[0mBase.PkgId, \u001b[90minput\u001b[39m::\u001b[0mString, \u001b[90mdepot_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mdl_load_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mload_path\u001b[39m::\u001b[0mVector\u001b[90m{String}\u001b[39m, \u001b[90mconcrete_deps\u001b[39m::\u001b[0mVector\u001b[90m{Pair{Base.PkgId, UInt64}}\u001b[39m, \u001b[90msource\u001b[39m::\u001b[0mNothing\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mloading.jl:1235\u001b[0m\n",
" [9] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
" [10] \u001b[0m\u001b[1meval\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90m./\u001b[39m\u001b[90;4mboot.jl:360\u001b[0m\u001b[90m [inlined]\u001b[39m\n",
" [11] \u001b[0m\u001b[1meval\u001b[22m\u001b[0m\u001b[1m(\u001b[22m\u001b[90mx\u001b[39m::\u001b[0mExpr\u001b[0m\u001b[1m)\u001b[22m\n",
"\u001b[90m @ \u001b[39m\u001b[90mBase.MainInclude\u001b[39m \u001b[90m./\u001b[39m\u001b[90;4mclient.jl:446\u001b[0m\n",
" [12] top-level scope\n",
"\u001b[90m @ \u001b[39m\u001b[90;4mnone:1\u001b[0m\n",
"in expression starting at /Users/gergaud/.julia/packages/DifferentialEquations/7WPQg/src/DifferentialEquations.jl:1\n"
]
},
{
"ename": "LoadError",
"evalue": "Failed to precompile DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa] to /Users/gergaud/.julia/compiled/v1.6/DifferentialEquations/jl_g1kMjp.",
"output_type": "error",
"traceback": [
"Failed to precompile DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa] to /Users/gergaud/.julia/compiled/v1.6/DifferentialEquations/jl_g1kMjp.",
"",
"Stacktrace:",
" [1] error(s::String)",
" @ Base ./error.jl:33",
" [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IJulia.IJuliaStdio{Base.PipeEndpoint}, internal_stdout::IJulia.IJuliaStdio{Base.PipeEndpoint}, ignore_loaded_modules::Bool)",
" @ Base ./loading.jl:1385",
" [3] compilecache(pkg::Base.PkgId, path::String)",
" @ Base ./loading.jl:1329",
" [4] _require(pkg::Base.PkgId)",
" @ Base ./loading.jl:1043",
" [5] require(uuidkey::Base.PkgId)",
" @ Base ./loading.jl:936",
" [6] require(into::Module, mod::Symbol)",
" @ Base ./loading.jl:923",
" [7] eval",
" @ ./boot.jl:360 [inlined]",
" [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
" @ Base ./loading.jl:1116"
]
}
],
"source": [
"# import Pkg; Pkg.add(\"DifferentialEquations\")\n",
"# Pkg.add(\"Plots\")\n",
"# import Pkg; Pkg.add(\"ODEInterfaceDiffEq\")\n",
"using DifferentialEquations\n",
"using Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pendule simple\n",
"### Introduction\n",
"On s'intéresse ici au pendule simple. Les principes physiques de la mécanique classique donnent comme équation qui régit l'évolution du mouvement\n",
"$$ ml^2\\ddot{\\alpha}(t)+mlg\\sin(\\alpha(t)) + k\\dot{\\alpha}(t)=0,$$\n",
"où $\\ddot{\\alpha}(t)$ désigne la dérivée seconde de l'angle $\\alpha$ par rapport au temps $t$. \n",
"\n",
"On prend ici comme variable d'état qui décrit le système $x(t)=(x_1(t),x_2(t))=(\\alpha(t), \\dot{\\alpha}(t))$. Le système différentiel du premier ordre que l'on obtient s'écrit alors\n",
"$$\n",
"\\left\\{\\begin{array}{l}\n",
"\\dot{x}_1(t) = x_2(t)\\\\\n",
"\\dot{x}_2(t) = -\\frac{g}{l}\\sin(x_1(t))-kx_2(t)\\\\\n",
"x_1(0) = x_{0,1}=\\alpha_0\\\\\n",
"x_2(0) = x_{0,2}=\\dot{\\alpha}_0\n",
"\\end{array}\\right.\n",
"$$\n",
"### Cas où la fonction second membre renvoie xpoint"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pendule (generic function with 1 method)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function pendule(x,p,t)\n",
"# second member of the IVP\n",
"# x : state\n",
"# real(2)\n",
"# p : parameter vector\n",
"# t : time, variable not use here\n",
"# real\n",
"# Output\n",
"# xpoint : vector of velocity\n",
"# same as x\n",
" g = p[1]; l = p[2]; k = p[3]\n",
" xpoint = similar(x)\n",
" xpoint[1] = x[2]\n",
" xpoint[2] = -(g/l)*sin(x[1]) - k*x[2]\n",
" return xpoint\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Main\n",
"#\n",
"g = 9.81; l = 10; k = 0;\n",
"p = [g,l,k] # constantes\n",
"theta0 = pi/3\n",
"t0 = 0.\n",
"tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période\n",
"tspan = (t0,tf) # instant initial et terminal\n",
"x0 = [theta0,0] # état initial\n",
"prob = ODEProblem(pendule,x0,tspan,p) # défini le problème en Julia\n",
"sol = solve(prob) # réalise l'intégration numérique\n",
"plot(sol, label = [\"x_1(t)\" \"x_2(t)\"]) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cas où xpoint est le premier argument modifié de la fonction second membre"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pendule1! (generic function with 1 method)"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function pendule1!(xpoint,x,p,t)\n",
"# second member of the IVP\n",
"# x : state\n",
"# real(2)\n",
"# p : parameter vector\n",
"# t : time, variable not use here\n",
"# real\n",
"# Output\n",
"# xpoint : vector of velocity\n",
"# same as x\n",
" g = p[1]; l = p[2]; k = p[3]\n",
" xpoint = similar(x)\n",
" xpoint[1] = x[2]\n",
" xpoint[2] = -(g/l)*sin(x[1]) - k*x[2]\n",
" return nothing\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Main\n",
"#\n",
"g = 9.81; l = 10; k = 0;\n",
"p = [g,l,k] # constantes\n",
"theta0 = pi/3\n",
"t0 = 0.\n",
"tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période\n",
"tspan = (t0,tf) # instant initial et terminal\n",
"x0 = [theta0,0] # état initial\n",
"prob = ODEProblem(pendule1!,x0,tspan,p) # défini le problème en Julia\n",
"sol = solve(prob) # réalise l'intégration numérique\n",
"plot(sol, label = [\"x_1(t)\" \"x_2(t)\"]) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"L'instruction `xpoint = similar(x)` **crée une variable locale xpoint** et on a perdu le paramètre formel `xpoint`."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pendule2! (generic function with 1 method)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function pendule2!(xpoint,x,p,t)\n",
"# second member of the IVP\n",
"# x : state\n",
"# real(2)\n",
"# p : parameter vector\n",
"# t : time, variable not use here\n",
"# real\n",
"# Output\n",
"# xpoint : vector of velocity\n",
"# same as x\n",
" g = p[1]; l = p[2]; k = p[3]\n",
" xpoint[1] = x[2]\n",
" xpoint[2] = -(g/l)*sin(x[1]) - k*x[2]\n",
" return nothing\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Main\n",
"#\n",
"g = 9.81; l = 10; k = 0.;\n",
"p = [g,l,k] # constantes\n",
"theta0 = pi/3\n",
"t0 = 0.\n",
"tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période\n",
"tspan = (t0,tf) # instant initial et terminal\n",
"x0 = [theta0,0] # état initial\n",
"prob = ODEProblem(pendule2!,x0,tspan,p) # défini le problème en Julia\n",
"sol = solve(prob) # réalise l'intégration numérique\n",
"plot(sol, label = [\"x_1(t)\" \"x_2(t)\"]) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Diagrammes de phases"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Main\n",
"#\n",
"g = 9.81; l = 10; k = 0.; # si k = 0.15 donc on a un amortissement\n",
"p = [g l k]\n",
"plot()\n",
"for theta0 in 0:(2*pi)/10:2*pi\n",
" theta0_princ = theta0\n",
" tf = 3*pi*sqrt(l/g)*(1 + theta0_princ^2/16 + theta0_princ^4/3072) # 2*approximation of the period\n",
" tspan = (0.0,tf)\n",
" x0 = [theta0 0]\n",
" prob = ODEProblem(pendule,x0,tspan,p)\n",
" sol = solve(prob)\n",
" plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
"end\n",
"theta0 = pi-10*eps()\n",
"x0 = [theta0 0]\n",
"tf = 70 # problem for tf=50 1/4 of the period!\n",
"tspan = (0.0,tf)\n",
"prob = ODEProblem(pendule,x0,tspan,p)\n",
"sol = solve(prob)\n",
"plot!(sol,vars=(1,2), xlims = (-2*pi,4*pi), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
"\n",
"theta0 = pi+10*eps()\n",
"x0 = [theta0 0]\n",
"tf = 70\n",
"tspan = (0.0,tf)\n",
"prob = ODEProblem(pendule,x0,tspan,p)\n",
"sol = solve(prob)\n",
"plot!(sol,vars=(1,2), xlims = (-2*pi,4*pi), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
"\n",
"# circulation case\n",
"for thetapoint0 in 0:0.2:2 \n",
" tf = 10\n",
" tspan = (0.,tf)\n",
" x0 = [-pi thetapoint0] # thetapoint0 > 0 so theta increases from -pi to ...\n",
" prob = ODEProblem(pendule,x0,tspan,p)\n",
" sol = solve(prob)\n",
" plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
"end\n",
"for thetapoint0 in -2:0.2:0\n",
" tf = 10\n",
" tspan = (0.,tf)\n",
" x0 = [3*pi thetapoint0] # thetapoint0 < 0 so theta decreases from 3pi to ...\n",
" prob = ODEProblem(pendule,x0,tspan,p)\n",
" sol = solve(prob)\n",
" plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
"end\n",
"plot!([-pi 0 pi 2*pi 3*pi], [0 0 0 0 0], seriestype=:scatter)\n",
"plot!(xlims = (-pi,3*pi))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modèle de Van der Pol\n",
"### Exemple de cycle limite\n",
"L'équation différentielle considérée est l'équation de Van der Pol\n",
"$$(IVP)\\left\\{\\begin{array}{l}\n",
"\\dot{y}_1(t)=y_2(t)\\\\\n",
"\\dot{y}_2(t)=(1-y_1^2(t))y_2(t)-y_1(t)\\\\\n",
"y_1(0)=2.00861986087484313650940188\\\\\n",
"y_2(0)=0\n",
"\\end{array}\\right.\n",
"$$\n",
"$t_f=T=6.6632868593231301896996820305$\n",
"\n",
"\n",
"La solution de ce problème de Cauchy est périodique de période $T$.\n",
"\n",
"Résoudre et de visualiser la solution pour les points de départ :\n",
"- x0 = [2.00861986087484313650940188,0]\n",
"- [x01 , 0] pour x01 in -2:0.4:0\n",
"- [0 , x02] pour x02 in 2:1:4\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"vdp (generic function with 1 method)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function vdp(x,p,t)\n",
"# Van der Pol model\n",
"# second member of the IVP\n",
"# x : state\n",
"# real(2)\n",
"# p : parameter vector\n",
"# t : time, variable not use here\n",
"# real\n",
"# Output\n",
"# xpoint : vector of velocity\n",
"# same as x\n",
"# To complete\n",
"# xpoint = similar(x)\n",
"# xpoint[1] = x[2]\n",
"# xpoint[2] = (1-x[1]^2)*x[2] - x[1]\n",
" xpoint = [x[2] , (1-x[1]^2)*x[2] - x[1]]\n",
" return xpoint\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x0 = [2.0086198608748433, 0.0], p = [1], tspan = (0.0, 6.66328685932313)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Main\n",
"#\n",
"tf = 6.6632868593231301896996820305\n",
"tspan = (0.0, tf)\n",
"x0 = [2.00861986087484313650940188,0]\n",
"mu = 1\n",
"p = [mu]\n",
"t0 = 0.;\n",
"tspan = (t0,tf)\n",
"println(\"x0 = $x0, p = $p, tspan = $tspan\")\n",
"prob = ODEProblem(vdp,x0,tspan,p)\n",
"sol = solve(prob)\n",
"plot(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false)\n",
"for x01 in -2:0.4:0 #4.5:4.5\n",
" x0 = [x01,0]\n",
" prob = ODEProblem(vdp,x0,tspan,p)\n",
" sol = solve(prob)#,force_dtmin=true)# \n",
" plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) \n",
"end\n",
"for x02 in 2:1:4\n",
" x0 = [0.,x02]\n",
" prob = ODEProblem(vdp,x0,tspan,p)\n",
" sol = solve(prob)#solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)\n",
" plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false)\n",
"end\n",
"plot!([0], [0], seriestype=:scatter) # point visualisation\n",
"plot!(xlims = (-2.5,5))\n",
" \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La solution périodique est ici ce qu'on appelle un **cycle limite**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Modèle de Lorentz\n",
"## Chaos\n",
"$$(IVP)\\left\\{\\begin{array}{l}\n",
"\\dot{x}_1(t)=-\\sigma x_1(t)+\\sigma x_2(t)\\\\\n",
"\\dot{x}_2(t)=-x_1(t)x_3(t)+rx_1(t)-x_2(t)\\\\\n",
"\\dot{x}_3(t)=x_1(t)x_2(t)-bx_3(t)\\\\\n",
"x_1(0)=-8\\\\\n",
"x_2(0)=8\\\\\n",
"x_3(0)=r-1\n",
"\\end{array}\\right.\n",
"$$\n",
"avec $\\sigma=10, r=28, b=8/3$.\n",
"\n",
"Calculer et visualisé la solution de ce problème\n"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lorentz (generic function with 1 method)"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function lorentz(x,p,t)\n",
"# Lorentz model\n",
"# second member of the IVP\n",
"# x : state\n",
"# real(3)\n",
"# p : parameter vector\n",
"# t : time, variable not use here\n",
"# real\n",
"# Output\n",
"# xpoint : vector of velocity\n",
"# same as x\n",
"# To complete\n",
" sigma = p[1]; rho = p[2]; beta = p[3];\n",
" xpoint = similar(x)\n",
" xpoint[1] = sigma*(x[2]-x[1])\n",
" xpoint[2] = x[1]*(rho-x[3]) - x[2]\n",
" xpoint[3] = x[1]*x[2] - beta*x[3]\n",
" return xpoint \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"xpoint = [160.0, -16.0, -136.0]\n"
]
}
],
"source": [
"# test de la fonction lorent!\n",
"sigma = 10.; rho = 28.; beta = 8/3;\n",
"p = [sigma rho beta]\n",
"#x0 = [1.1,0.0,0.0]\n",
"x0 = [-8, 8, rho-1]\n",
"xpoint = x0\n",
"xpoint = lorentz(x0,p,0)\n",
"println(\"xpoint = $xpoint\")"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Main\n",
"#\n",
"x0 = [-8, 8, rho-1]\n",
"tspan = (0.0,100.0)\n",
"prob = ODEProblem(lorentz,x0,tspan,p)\n",
"sol = solve(prob)\n",
"plot(sol,vars=(1,2,3), xlabel = \"x_1\", ylabel = \"x_2\", zlabel = \"x_3\", legend = false)\n",
"plot!([x0[1]], [x0[2]], [x0[3]], seriestype=:scatter) # point visualisation\n",
"#annotate!([x0[1],x0[2],x0[3],text(\". Point de départ\", 10,:left)])"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p1 = plot(sol,vars=(0,1),ylabel = \"x_1(t)\", legend = false)\n",
"p2 = plot(sol,vars=(0,2),ylabel = \"x_2(t)\", legend = false)\n",
"p3 = plot(sol,vars=(0,3),xlabel = \"t\", ylabel = \"x_3(t)\", legend = false)\n",
"plot(p1,p2,p3,layout=(3,1))#,legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Attention à bien utiliser les codes\n",
"Nous allons ici résoudre le problème à valeur initiale\n",
"$$(IVP)\\left\\{\\begin{array}{l}\n",
"\\dot{x}_1(t)=x_1(t)+x_2(t)+\\sin t\\\\\n",
"\\dot{x}_2(t)=-x_1(t)+3x_2(t)\\\\\n",
"x_1(0)=-9/25\\\\\n",
"x_2(0)=-4/25,\n",
"\\end{array}\\right.\n",
"$$\n",
"dont la solution est\n",
"\\begin{align*}\n",
"x_1(t)&= (-1/25)(13\\sin t+9\\cos t)\\\\\n",
"x_2(t)&= (-1/25)(3\\sin t+4\\cos t).\n",
"\\end{align*}"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function exemple1(x,p,t)\n",
"# second member of the IVP\n",
"# Input\n",
"# x : state\n",
"# real(2)\n",
"# p : parameter vector\n",
"# t : time\n",
"# real\n",
"# Output\n",
"# xpoint : vector of velocity\n",
"# same as x\n",
" xpoint = similar(x) # xpoint est un objet de même type que x\n",
" xpoint[1] = x[1] + x[2] + sin(t)\n",
" xpoint[2] = -x[1] + 3*x[2] \n",
"# xpoint = [x[1]+x[2]+sin(t), -x[1] + 3*x[2]]\n",
" return xpoint \n",
"end\n",
"p = []\n",
"t0 = 0.; tf = 8\n",
"tspan = (t0,tf)\n",
"x0 = [-9/25 , -4/25]\n",
"prob = ODEProblem(exemple1,x0,tspan,p) # défini le problème en Julia\n",
"sol = solve(prob, DP5()) # réalise l'intégration numérique\n",
" # avec ode45 de matlab\n",
"p1 = plot(sol, label = [\"x_1(t)\" \"x_2(t)\"], title = \"ode45 de Matlab\") #, lw = 0.5) # lw = linewidth\n",
"sol = solve(prob, DP5(), reltol = 1.e-6, abstol = 1.e-9)\n",
"p2 = plot(sol, label = [\"x_1(t)\" \"x_2(t)\"], title = \"reltol=1.e-6, abstol=1.e-9\") #, lw = 0.5)\n",
"sol = solve(prob, DP5(), reltol = 1.e-10, abstol = 1.e-15)\n",
"p3 = plot(sol, label = [\"x_1(t)\" \"x_2(t)\"], title = \"reltol=1.e-10, abstol=1.e-16\") #, lw = 0.5)\n",
"T = t0:(tf-t0)/100:tf\n",
"sinT = map(sin,T) # opération vectorielle\n",
"cosT = map(cos,T)\n",
"p4 = plot(T,[(-1/25)*(13*sinT+9*cosT) (-1/25)*(3*sinT+4*cosT)], label = [\"x_1(t)\" \"x_2(t)\"], xlabel = \"t\", title = \"Solution exacte\")\n",
"plot(p1,p2,p3,p4)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.6.3",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}