{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Time-dependent PDEs\n",
"\n",
"Copyright (C) 2010-2020 Luke Olson
\n",
"Copyright (C) 2020 Andreas Kloeckner\n",
"\n",
"\n",
"MIT License
\n",
"Permission is hereby granted, free of charge, to any person obtaining a copy\n",
"of this software and associated documentation files (the \"Software\"), to deal\n",
"in the Software without restriction, including without limitation the rights\n",
"to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
"copies of the Software, and to permit persons to whom the Software is\n",
"furnished to do so, subject to the following conditions:\n",
"\n",
"The above copyright notice and this permission notice shall be included in\n",
"all copies or substantial portions of the Software.\n",
"\n",
"THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
"IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
"FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
"AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
"LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
"OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
"THE SOFTWARE.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as pt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"def rk4_step(y, t, h, f):\n",
" k1 = f(t, y)\n",
" k2 = f(t+h/2, y + h/2*k1)\n",
" k3 = f(t+h/2, y + h/2*k2)\n",
" k4 = f(t+h, y + h*k3)\n",
" return y + h/6*(k1 + 2*k2 + 2*k3 + k4)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"mesh = np.linspace(0, 1, 200)\n",
"dx = mesh[1]-mesh[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(all of the PDEs below use *periodic* boundary conditions: $u(0)=u(1)$)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Advection equation:** $u_t+u_x=0$\n",
"\n",
"Equivalent: $u_t=-u_x$"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"def f_advection(t, u):\n",
" du = (np.roll(u, -1, axis=-1) - np.roll(u, 1, axis=-1))/(2*dx)\n",
" return -du"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Heat equation:** $u_t=u_{xx}$"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"def f_heat(t, u):\n",
" d2u = (\n",
" np.roll(u, -1, axis=-1)\n",
" - 2*u\n",
" + np.roll(u, 1, axis=-1))/(dx**2)\n",
" return d2u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Wave equation:** $u_{tt}=u_{xx}$\n",
"\n",
"NOTE: Two time derivatives $\\rightarrow$ convert to first order ODE.\n",
"\n",
"$$u_t=v$$\n",
"$$v_t=u_{xx}$$\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"def f_wave(t, w):\n",
" u, v = w\n",
" d2u = (\n",
" np.roll(u, -1, axis=-1)\n",
" - 2*u\n",
" + np.roll(u, 1, axis=-1))/(dx**2)\n",
" return np.array([v, d2u])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Initial condition**"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"data": {
"text/plain": [
"(1, 200)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"current_t = 0\n",
"\n",
"#current_u = np.sin(2*np.pi*mesh)*0.5+1\n",
"#current_u = (mesh > 0.3) & (mesh < 0.7)\n",
"#current_u = (mesh > 0.45) & (mesh < 0.55)\n",
"current_u = np.exp(-(mesh-0.5)**2*150)\n",
"#current_u = 2*np.abs(mesh-0.5)\n",
"\n",
"current_u = np.array([current_u], dtype=np.float64)\n",
"current_u.shape"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"data": {
"text/plain": [
"(2, 200)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Add a second component if needed (for wave equation)\n",
"current_u = np.vstack([current_u,np.zeros(len(mesh))])\n",
"current_u.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Time loop**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run this cell many times in place (using Ctrl-Enter):"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD7CAYAAACMlyg3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVXX9x/HXR3AtYTQXFJEJwSW3AdO0CEZQwS2s3DDR\nKbcWUyoNKg1xS8USjKyfqY1LimuapoIpFzBU3MYFAU1CQcxcMM0VmM/vj+8dzjgOw8zcc+855973\n8/G4D+fMPXPOxw/33s89388532PujoiIVKa1kg5ARESSoyIgIlLBVARERCqYioCISAVTERARqWAq\nAiIiFaxr0gE0Z2Y6X1VEpBPc3Trzd6k7EnB3PdwZN25c4jGk5aFcKBfKRduPQqSuCEiwaNGipENI\nDeUiolxElIt4xFIEzOxKM3vNzJ5ezfNHmdlTZtZgZg+a2c5x7FdERAoT15HAn4BhbTy/EBjk7jXA\nucAfY9pv2aqrq0s6hNRQLiLKRUS5iIcVOp60akNmvYE73X2XNaxXBTzj7r1aec7jikdEpFKYGZ6h\nxvDxwD0J7DdTcrlc0iGkhnIRUS4iykU8SnqKqJntDXwbGLi6derq6qiurgagqqqKmpoaamtrgegf\nXcuVtdwkLfEkudzQ0JCqeJJcbmhoSFU8pVzO5XLU19cDrPq87KySDQeZ2S7ArcBwd39xNetoOEhE\npIPSMhxk+cennzDbmlAARq2uAIiISOnFdYro9cBsYFsze9nMvm1mJ5nZiflVzgQ2Bi4zsyfNbE4c\n+y1nLYdCKplyEVEuIspFPGLpCbj7UWt4/gTghDj2JSIi8YmtJxAH9QRERDouLT0BERHJGBWBlNJ4\nZ0S5iCgXEeUiHioCIiIVTD0BEZGMU09AREQ6RUUgpTTeGVEuIspFRLmIh4qAiEgFU09ARCTj1BMQ\nEZFOURFIKY13RpSLiHIRUS7ioSIgIlLB1BMQEck49QRERKRTVARSSuOdEeUiolxElIt4qAiIiFQw\n9QRERDJOPQEREekUFYGU0nhnRLmIKBcR5SIeKgIiIhUslp6AmV0JHAS85u67rGadS4H9gfeAOndv\naGUd9QRERDooDT2BPwHDVvekme0PbOPu/YCTgD/EtF8RESlALEXA3R8ElrWxygjgmvy6jwDdzWzz\nOPZdrjTeGVEuIspFRLmIR6l6Aj2Bxc2WX8n/TkREEhTbdQJm1hu4s7WegJndBZzv7rPzy38HTnf3\nJ1usp56AiEgHFdIT6Bp3MKuxBOjVbHkrYGlrK9bV1VFdXQ1AVVUVNTU11NbWAtHhn5a1rGUtV/Jy\nLpejvr4eYNXnZWfFeSRQTTgS2LmV5w4AfuDuB5rZnsBEd9+zlfV0JJCXy+VW/eNXOuUiolxElItI\n4kcCZnY9UAt8zsxeBsYB6wDu7pe7+91mdoCZ/ZNwiui349iviIgURnMHiYhkXBquExARkQxSEUip\npiaQKBfNKRcR5SIeKgIiIhVMPQERkYxTT0BERDpFRSClNN4ZUS4iykVEuYiHioCISAVTT0BEJOPU\nExARkU5REUgpjXdGlIuIchFRLuKhIiAiUsHUExARyTj1BEREpFNUBFJK450R5SKiXESUi3ioCIiI\nVDD1BEREMk49ARER6RQVgZTSeGdEuYgoFxHlIh4qAiIiFUw9ARGRjFNPQEREOiWWImBmw81svpk9\nb2ZjWnm+l5k9YGZPmFmDme0fx37LmcY7I8pFRLmIKBfxKLgImNlawGRgGLAjMNLMtm+x2hnAje4+\nABgJXFbofkVEpHAF9wTMbE9gnLvvn18eC7i7X9hsnd8DC919gpntBUxw94GtbEs9ARGRDiqkJ9A1\nhv33BBY3W14C7NFinfHANDM7BdgA2CeG/YqISIHi6Am0Vn1afp0fCfzJ3XsBBwLXxbDfsqbxzohy\nEVEuIspFPOI4ElgCbN1seStgaYt1jiP0DHD3h81sPTPbxN3faLmxuro6qqurAaiqqqKmpoba2log\n+kfXcmUtN0lLPEkuNzQ0pCqeJJcbGhpSFU8pl3O5HPX19QCrPi87K46eQBdgATAUeBWYA4x093nN\n1vkbcJO7X21mOwD3uftWrWxLPQERkQ5K9DoBd18JnAxMA+YCU9x9npmNN7OD8qudBpxgZg3An4Fj\nC92viIgULpbrBNz9Xnffzt37ufsF+d+Nc/e78j/Pc/eB7l7j7gPc/f449lvOWg6FVDLlIqJcRJSL\neOiKYRGRCqa5g0REMk5zB4mISKeoCKSUxjsjykVEuYgoF/FQERARqWDqCYiIZJx6AiIi0ikqAiml\n8c6IchFRLiLKRTxUBEREKph6AiIiGaeegIiIdIqKQEppvDOiXESUi4hyEQ8VARGRCqaegIhIxqkn\nICIinaIikFIa74woFxHlIqJcxENFQESkgqknICKSceoJiIhIp6gIpJTGOyPKRUS5iCgX8VAREBGp\nYLH0BMxsODCRUFSudPcLW1nncGAc0Ag85e5Ht7KOegIiIh1USE+g4CJgZmsBzwNDgaXAo8CR7j6/\n2Tp9gRuBvd39HTPbxN3faGVbKgIiIh2UdGN4D+AFd3/J3ZcDU4ARLdY5Afidu78D0FoBkE/SeGdE\nuYgoFxHlIh5xFIGewOJmy0vyv2tuW2A7M3vQzGab2bAY9isiIgXqGsM2WjsEaTmm0xXoCwwCtgZm\nmdmOTUcGzdXV1VFdXQ1AVVUVNTU11NbWAlHlr4Tl2traVMWj5fQsN0lLPEktN/0uLfGUcjmXy1Ff\nXw+w6vOys+LoCewJnOXuw/PLYwFv3hw2s98DD7n7NfnlvwNj3P3xFttST0BEpIOS7gk8CvQ1s95m\ntg5wJPDXFuvcDgwBMLNNgH7Awhj2XbZafuurZMpFRLmIKBfxKLgIuPtK4GRgGjAXmOLu88xsvJkd\nlF9nKvCmmc0F7gdOc/dlhe5bREQKo7mDREQyLunhIBERySgVgZTSeGdEuYgoFxHlIh4qAiIiFUw9\nARGRjFNPQEREOkVFIKU03hlRLiLKRUS5iIeKgIhIBVNPQEQk49QTEBGRTlERSKlCxjuXLIHZs+Gh\nh+Cjj+KLKSka+42UQy5efx1mzYKHH4b//a/z2ymHXKSBikAZuf9+GDgQdt0VTjsNTj4ZNt8cjj8+\nvPFEkjRnDgwZAv36wdix8MMfQs+e8K1vweLFa/57KQ4VgZRqPmf6mixfHj7w6+rg1FPh1VfDkcDj\nj8OCBVBVBTvvDHffXbRwi6ojuSh3WcyFO/ziF/C1r8Gxx8K//w3/+Ac8+ii8+CJstx0MGAA33NCx\n7WYxF2mkxnDGffwxjBwJH3wA118fPvBbM3s2fP3rMHkyHHZYaWOUyrVyJXz3uzB3LtxxB2y6aevr\nPfUUHHwwnHEGnHhiaWMsB2oMl6H2jHe6wzHHhCOBv/xl9QUA4Mtfhvvug1NOgXvvjS/OUtDYbyRr\nufjRj+CFF2DatNUXAAhDmNOnw3nnwbXXtm/bWctFWsVxe0lJyMUXh8PpWbNg3XXXvP4uu8BNN8Gh\nh4bD8b59ix+jVK76+vCFY84c+Oxn17z+NtvA3/4Ge+8NO+0E/fsXPURBw0GZNWMGHHFEeINtvXXH\n/vayy+D//i+Mya6zTnHik8r27LPhw3zGDPjCFzr2tzfdBGPGQEMDdO9enPjKTSHDQSoCGfT++6HR\ne8klodnWUe7h7wYMgPHj449PKtuKFbDXXmFs/4QTOreNE08Es/BlRdZMPYEy1NZ455lnwp57dq4A\nQPTm+v3vw7ettNPYbyQLufj1r0N/6vjjO7+NCRPC2WzTp69+nSzkIgvUE8iYp56C664LZ1sUYsst\n4fzzw6mls2aFwiBSqMWL4aKL4LHHCntNde8ehi1POikMLWnYsng0HJQh7jB0aDjF83vfK3x7K1fC\n7rvDT38KRx5Z+PZEvvUt6NMHzjknnu3tvz8MGwajR8ezvXKlnkCFuOOOcNFNQwN0jekYbuZMGDUK\n5s+H9dePZ5tSmR56KHxBmT+/fWcDtcdzz8HgwWGbn/tcPNssR4n3BMxsuJnNN7PnzWxMG+sdamaN\nZjYgjv2Ws5bjnStXhkvtL744vgIAMGgQ7LZbOPROK439RtKaC/fw+jznnPgKAIQziw47LFw/0FJa\nc5E1BRcBM1sLmAwMA3YERprZ9q2s91ngh8DDhe6zEt1yC2y0UTg0jtvZZ4dGXCGTeUllu//+MB3E\nqFHxb/vMM8M1B0uXxr9tiWE4yMz2BMa5+/755bGAu/uFLda7BLgPOB34ibs/0cq2NBzUisbGcKHX\nxRfD8OHF2cdRR4ULdH7+8+JsX8qXezgldPTo4vWWTjsNPvwwTHsin5b0cFBPoPkcgEvyv1vFzGqA\nrdw9o1OYJeu222CDDYpzFNBk3Lhw3YGOBqSjpk0Lr5vDDy/ePsaMCRPMLVlSvH1UqjhGl1urPqu+\nzpuZAZcAx67hbwCoq6ujuroagKqqKmpqalbNFtg0BlgJy00/NzbCOefUcv75MGNG8fa33XbwhS/k\nGDsWJk9O/v+/+XLLnCQdT5LLDQ0NjM6fKpOGeADOP7+Wn/8cZs4s3v423RSGDMnxk5/AjTeG5ydO\nnFjRnw/19fUAqz4vO83dC3oAewL3NlseC4xpttwN+A+wEPgX8AHhaGFAK9tyCaZPn+7u7rfd5j5g\ngHtjY/H3OWeOe69e7h9/XPx9dURTLiR9uZg1y71PH/fly4u/r5dfdt94Y/e33grLactFkvKfnZ36\nDI+jJ9AFWAAMBV4F5gAj3X3eatafDvzY3Z9s5TkvNJ5y4h6mdjjrLBgxojT7HDIEvv3t4jT4pPwc\neGB4bZZq+udjjw33H1Dv6pMS7Qm4+0rgZGAaMBeY4u7zzGy8mR3U2p/QxnCQRP7+9zAPS2enh+iM\nMWPCFZ+qxbIm8+eHK4OPOaZ0+zz9dLj00tAklnjEcp2Au9/r7tu5ez93vyD/u3Huflcr6w7xVs4M\nkk/K5XJMnBjmYy/llA777QddusA995Run2vSvDdQ6dKUi0svDTeMWW+90u1zp53gi1+Eq69OVy6y\nTBPIpdTLL4dvWUcdVdr9moVpJC68cM3rSuVatiycrRPH9CUd9dOfhtOlV64s/b7LkYpASj3ySC0n\nnVTab1lNDj88FKFHHin9vlvTdHaEpCcXV1wRbgfZo0fp9/3Vr8Imm8CyZbWl33kZ0txBKbRsWbjL\n0ty5sMUWycRwySXhhjUdvfm3lL8VK8Lr87bbwpQjSbj11jBl9ezZyew/bZK+WExidsUVsNtuucQK\nAMB3vgNTp6bj4hyN/UbSkIvbb4devZIrAACHHAILF+ZSc7SaZSoCKbNiRbg0/tBDk42je/dw1ocu\n05eWJk5MfmrnLl3gm98MR6xSGA0Hpcwtt8CkSeFGL0l78cVwB7NFi+Azn0k6GkmDxx+Hb3wjvDbi\nnM22M955Bz7/eXjyyY7fZ7vcaDiojKThW1aTbbaBgQPhmmuSjkTSYtKkcDe6pAsAQLdu4eKx3/42\n6UiyTUUgRR59NIzBjxiRjrFfCAVp0qQwh1FS0pKLNEgyF6++CnfeWdi9g+OUy+U45RS46ipNfFgI\nFYEUmTQJfvjDdHzLajJoUJjBdOrUpCORpP3hDzByZLivRVpUV8Pee8Of/pR0JNmlnkBKLF0aroZc\nuBCqqpKO5pOuuSbc3H7atKQjkaR8+CH07g0zZsD2n7plVLL+8Y8wLLRgQWgYVyL1BMrAZZeFm3Sn\nrQAAHHEEPPNMuG5BKtOUKdC/f/oKAMCXvwwbbwx3fWqSGmkPFYEU+OADuPxyOOWU6HdpGgdfd134\n/vdD0zoJacpF0pLIhXu6Tlho0pQLM/jxj+E3v0k2nqxSEUiBP/8ZvvQl6Ncv6UhW76STwumrr7+e\ndCRSajNnhuGg/fZLOpLV++Y3w1DqE5qassPUE0iYO+y8c2gKDx2adDRtO/740Ig744ykI5FS+vrX\nQwFIYrK4jrjoojBsee21SUdSeoX0BFQEEvb3v4fpop9+urRTRnfGM8+E+xwvWgTrrJN0NFIK//oX\n7L47vPRS+i8YXLYM+vQJvastt0w6mtJSYzjDmsZaWxaANI6D77wz7Lgj3HRTafebxlwkpdS5mDw5\nzCOVxgLQMhcbbRROrvjd75KJJ6tUBBL0/PPhArFS3zOgEKNHh/laKuyArSK9+y7U18MPfpB0JO13\n6qnwxz/C++8nHUl2aDgoQd//fji17dxzk46k/RobYYcdwkynX/1q0tFIMU2aBA8+CDffnHQkHTNi\nBBxwQDiZoVKoJ5BBr78ebpg9bx5svnnS0XTMZZeFXsZttyUdiRTLihXQty/ceGM4cy1Lcrlw28vn\nnoO1KmSsQz2BDJo8GQ47bPUFIM3j4MccE04bXLiwNPtLcy5KrVS5uOWWMDNnmgvA6nIxeDCsvz7c\ne29p48mqWIqAmQ03s/lm9ryZjWnl+R+Z2VwzazCz+8ysVxz7zar33oPf/x5+8pOkI+mcz34WjjtO\nszeWK3eYMAFOPz3pSDrHLJxxp3sNtE/Bw0FmthbwPDAUWAo8Chzp7vObrTMYeMTdPzSz7wK17n5k\nK9uqiOGgyZPhgQeyPZyyeDHU1IRTCLt1SzoaidP06aFfNXdudodTPv44XNMydWo4q63cJT0ctAfw\ngru/5O7LgSnAiOYruPsMd/8wv/gw0DOG/WbSihXh8vaf/jTpSArTqxfsu69mbyxHEyaEo9SsFgAI\n17H84AfJTXWSJXH8M/cEFjdbXkLbH/LHAffEsN9MuvVW2GqrcMeutmRhHLzpXgMrVxZ3P1nIRakU\nOxfPPhvu1HX00UXdTSzWlIuTTgpH2//5T2niyao4ikBrhyCtjumY2dHAbsCEGPabOe5w4YXZHWtt\nac89Q2P7L39JOhKJy4QJ4c5h662XdCSF22QTOPxw3Sd7TeK4fckSoPkdPrci9AY+wcz2AX4GDMoP\nG7Wqrq6O6upqAKqqqqipqaG2thaIKn9Wl887L8c778CBB655/dra2sTjbc/y174G48fX8o1vwMyZ\nycdTCctN4t7+tdfmuP12ePnldP3/rm656XdtrT9oEJx6ai2jR8PTT6cr/kKWc7kc9fX1AKs+Lzsr\njsZwF2ABoTH8KjAHGOnu85qt0x+4GRjm7i+2sa2ybQw3Nob52M89Fw4+OOlo4uMejghOOy2c8irZ\nNWpUuHal3CYIPOEE2GwzOO+8pCMpnkQbw+6+EjgZmAbMBaa4+zwzG29mB+VXuwj4DHCzmT1pZrcX\nut+sufXWMC//QQeteV3Izji4GZx1VngUqzeQlVyUQrFyMX9+OJOm+T0t0q69uTjjjHBrzDfeKG48\nWRVL/9/d73X37dy9n7tfkP/dOHe/K//zvu6+hbsPcPf+7n5IHPvNipUrYdw4OPvs9M8U2hnDh4fT\nRLM2vYBEzj47NPrL8XTf3r3D3fEuuijpSNJJ00aUwJ//HKZaePDB8iwCEO4/fOqp4eySSr3Pa1bN\nnRtu1v7ii7DhhklHUxxLlsAuu4SpJHr0SDqa+CV9nYC0YcUKGD++fI8Cmuy7b5gMb8qUpCORjjr7\n7HBdQLkWAAinZY8aFc7Ok09SESiyq66Cnj1hyJCO/V3WxsHNQtP7zDPho4/i3XbWclFMcefiscfC\nPFBZmi66SUdz8bOfwTXXhKvcJaIiUET//S/88pfw61+X91FAk733Djed0ZxC2eAebtB+9tlhPqhy\n16NH6HuM+dTsZpVNPYEiGjMmTBl91VVJR1I6CxbAV74SpsjedNOko5G23HZbOKvryScrp4/z/vuw\n/fZw/fUwcGDS0cRH9xNIofnzw4vsmWdgiy2Sjqa0Tj01vNn++MekI5HVef/9cNR2xRUwdGjS0ZTW\nDTeE3sBjj0HXOC6XTQE1hlPGPdzU4pe/7HwByPI4+Nlnw913w+zZ8Wwvy7mIW1y5OPfccJFflgtA\nZ3Nx5JHhKPXSS+ONJ6tUBIrgmmvgf//LZrMtDt27hz7I974Hy1c7QYgk5bnnwlHab36TdCTJMAun\nbJ9/Prz8ctLRJE/DQTF75ZUwPcTUqeG/lcodDjww3Jlq3Liko5EmK1bAl78M3/lOOFqtZL/6Vbiv\nx9Sp2Z42G9QTSA132H//0Bg988yko0leU0G8914YMCDpaATCt9/p08PFfZVwxlpbVqwIfbtRo7J/\n1K6eQEr89rfw5pswdmzh2yqHcfCePcOQw6hR4ZaanVUOuYhLIbl47LFwk5UrryyPAlDo66JrV7j6\n6nCG1LPPxhJSJqkIxOSRR0Kz7cYbYe21k44mPb71LfjiF7P/TSvr3n47zK1/2WXhBvISbLcdXHxx\nmAH3f/9LOppkaDgoBq+9Fsa+L7kEvv71pKNJn/feg913Dzf/PuGEpKOpPI2N4XW59da6kG91jjsu\nFIEbbshmf0A9gQR98EG4UnbYsDBHkLTu+efhq18Nb7KOTqEhhRkzBh5+OPQB1l036WjS6cMPw+ty\nn33CKc5Zo55AQlasCMMdffqEccU4lds4+LbbhgIwcmTHx1/LLReF6Ggu/vCHcC+LpvtZlJM4Xxfr\nrQe33w7XXRcuoKskZXK9XOmtXAnHHhuOBG64oTwabcU2ZEhoTA4bFs5Q2XbbpCMqb9deG+6mlcuF\n++1K2zbbLJwuuvfeoSgcfXTSEZWGhoM64aOPwhkvb70Fd94J66+fdETZctVV4Wrqe+6BnXdOOpry\ndMUVIcf33w877JB0NNny3HNhWOiss+DEE5OOpn0KGQ7SkUAHvfVWOJNgo43grrvCNwbpmO98JxTO\noUPhppug2X3DpUDu4dv/lVeGIwAdbXXcF74AM2aEO+a98kq42DGLzeL2KuP/tfg98UQ43bF//3Aq\naDELQLmPg48cGYbRjjgiXEvQ1gFgueeiI9rKxTvvwKGHhqPTf/yj/AtAMV8X/fqFHN5/P4wYAcuW\nFW1XiVMRaIfly+Gcc8I3gwsuCOcVV8rUu8U0dGi4vuKGG8KV1kuWJB1Rdj3wQLh94mabhZvEbLll\n0hFlX48eIa99+4Zhy7vvTjqi4lBPoA3uoVH04x9DdTVcfnm4TZ3Ea/nyUFwnTYLTTgs3/tAwW/ss\nWRKuUM/lwuvzgAOSjqg8PfAAHH881NTAhAmwzTZJR/RJiZ8iambDzWy+mT1vZp+6b4+ZrWNmU8zs\nBTN7yMxSfc3i8uVwxx1hDqDRo8Pc43/7mwpAsay9dphr6aGHwvnsffuGs4j++9+kI0uvhQvDVdi7\n7BK+oMyfrwJQTEOGhIbx7ruHC0Pr6qChIemo4lFwETCztYDJwDBgR2CkmW3fYrXjgLfcvR8wEbio\n0P3GzT2M+Y8eHT7sL7oo/Dx3Lhx8cOlPAa3EcfB+/cK52nfcEQpC796hZ/CrX+U0JTVhXHrs2Bz7\n7Rc+iDbcMNzB7dxzK+P2kC2V+j2y3nrhPsX//Gd4rR5yCOy6a5g2/ZVXShpKrOI4O2gP4AV3fwnA\nzKYAI4D5zdYZATRNKHwLoWgk6t13w7enOXPCofSMGdCtGxx1VGgI9e2bdISVa7fdQuP9rbfg5pvD\nVAcXXAB77QWDBoWZH3faCTbeOOlIi6exERYvhqeeCq/NGTPCrTt33TUcAfz1rxoyS0pVFfziF6Eg\nzJwZ7h+y006hh1BbC4MHh2GjbbbJxjxiBfcEzOybwDB3PzG/fDSwh7uf0mydZ/LrLM0vvwB8yd3f\narGtdvcE3MOjsRE+/jjcLq/l4733wj1+X3sN/v3v8N9XXw1TGLz+epg8qn//cHHI4MGaWCvN3nwT\nZs0Kb7rZs8Oh+frrh/vF9u4dGqFbbBEao926hW/J3bqFxwYbhDdj80eXLsU9unMPFxSuWPHJx8cf\nhzlq3nknevz3v+H1+eqrsHQpvPhi+MDv3j18uAwaFF6fu+9eflf9louVK+Hpp8MXypkzw1XxixeH\nQrDNNmFG3S23DIWiqip6bXbvHo7i1lmnsNdoonMHmdmhwH4tisDu7n5qs3Weza/TVAT+mV9nWYtt\n+ec+5zQ2ssaHe0iQWXhjbLBBeKy/fvTzBhuE28j16AGbbx7+26NHOJTr3Vtn+GSZe/jAnD8/vNmW\nLg0fov/5TzjKe/fd6EP2/fdDn6fpsWJFeNN27Rq94drzZmtrndY+9Lt0Cfto/lh77fCmb/oQaCpW\nPXpEhaxPn1DcunePL19Seh98EIr5okXh9fnKK6HYNxX+pv++++4nX59Nj8bGTxeE1T3eeCPZi8WW\nAM2/Q28FLG2xzmKgF7DUzLoA3VoWgCb77FNH797VmEH37lXU1NQweHAta60Fs2blWGst2HvvWsxg\nxowcALX5q42axgjXtNynT8fWT2K5+XhnGuJJcrnpd82fN4MXXsjRpQvU1XV8+42NcP/9OVasgL32\nCs8/+GB4fuDATy+7t/08wEMPhXiGDKmlSxeYObNz/79f+tLqn29oaGD06NEd/v8tx+WJEydSU1OT\nmnhaLj/ySFg+5JDO/f0DD+Tyd4KrZeXK8PnnDl/5Si2zZuW44YZ6AHr1qubii+m0OI4EugALgKHA\nq8AcYKS7z2u2zveBndz9+2Z2JHCIux/ZyrZSdYpoknK53KoXQ6VTLiLKRUS5iCQ+lbSZDQcmEc42\nutLdLzCz8cCj7n6Xma0LXAv0B94EjnT3Ra1sR0VARKSDEi8CcVEREBHpuMQvFpP4NR8Pr3TKRUS5\niCgX8VAREBGpYBoOEhHJOA0HiYhIp6gIpJTGOyPKRUS5iCgX8VAREBGpYOoJiIhknHoCIiLSKSoC\nKaXxzohyEVEuIspFPFQEREQqmHoCIiIZp56AiIh0iopASmm8M6JcRJSLiHIRDxUBEZEKpp6AiEjG\nqScgIiKdoiKQUhrvjCgXEeUiolzEQ0VARKSCqScgIpJx6gmIiEinqAiklMY7I8pFRLmIKBfxKKgI\nmNlGZjbNzBaY2VQz697KOrua2Wwze8bMGszs8EL2KSIi8SmoJ2BmFwJvuvtFZjYG2Mjdx7ZYpy/g\n7v6imW0BPA5s7+7vtLI99QRERDqokJ5AoUVgPjDY3V8zsx5Azt23X8PfNADfdPcXW3lORUBEpIOS\nbAxv5u4kfnkcAAAFGElEQVSvAbj7v4FN21rZzPYA1m6tAMgnabwzolxElIuIchGPrmtawczuAzZv\n/ivAgTM6sqP8UNA1wKi21qurq6O6uhqAqqoqampqqK2tBaJ/dC1X1nKTtMST5HJDQ0Oq4klyuaGh\nIVXxlHI5l8tRX18PsOrzsrMKHQ6aB9Q2Gw6a7u47tLLehkAOOM/db2tjexoOEhHpoCSHg/4K1OV/\nPha4o+UKZrY2cDtwdVsFQERESq/QInAhsK+ZLQD2AS4AMLPdzOzy/DqHAwOBOjN70syeMLNdCtxv\n2Ws5FFLJlIuIchFRLuJRUBFw97fcfR93387d93X3t/O/f9zdT8z//Gd3X9fdB7h7//x/n44j+HLW\nNN4pykVzykVEuYiHrhhOqbfffjvpEFJDuYgoFxHlIh4qAiIiFUxFIKUWLVqUdAipoVxElIuIchGP\n1E0lnXQMIiJZlMi0ESIikm0aDhIRqWAqAiIiFSyRImBmw81svpk9n5+CuuXz65jZFDN7wcweMrOt\nk4izFNqRix+Z2dz8vRjuM7NeScRZCmvKRbP1DjWzRjMbUMr4Sqk9uTCzw/OvjWfM7LpSx1gq7XiP\n9DKzB/IXojaY2f5JxFlsZnalmb1mZqu9zsrMLs1/bjaYWU27NuzuJX0QCs8/gd7A2kAD4f4Czdf5\nHnBZ/ucjgCmljjNFuRgMrJf/+buVnIv8ep8FZgCzgQFJx53g66Iv4d4c3fLLmyQdd4K5+D/gpPzP\nOwD/SjruIuViIFADPL2a5/cH/pb/+UvAw+3ZbhJHAnsAL7j7S+6+HJgCjGixzgjg6vzPtwBDSxhf\nKa0xF+4+w90/zC8+DPQscYyl0p7XBcA5hOlKPiplcCXWnlycAPzO8zdncvc3ShxjqbQnF41At/zP\nVcArJYyvZNz9QWBZG6uMIMzUjLs/AnQ3s83bWB9IZjioJ7C42fISPv3Btmodd18JvG1mG5cmvJJq\nTy6aOw64p6gRJWeNucgf3m7l7neXMrAEtOd1sS2wnZk9mL9967CSRVda7cnFeGCUmS0G7gJ+WKLY\n0qZlrl6hHV8a13g/gSJo7VzWlueptlzHWlmnHLQnF2FFs6OB3QjDQ+WozVyYmQGXEGarbetvykF7\nXhddCUNCg4CtgVlmtqO3ctvWjGtPLkYCf3L3S8xsT+A6YMeiR5Y+7f48aS6JI4ElhBdtk62ApS3W\nWQz0AjCzLoRxz7YOg7KqPbnAzPYBfgYcnD8kLkdrysWGhDd2zsz+BewJ3FGmzeH2vC6WAHe4e6O7\nLwIWAP1KE15JtScXxwE3Abj7w8B6ZrZJacJLlSXkPzfzWv08aSmJIvAo0NfMepvZOsCRhPsSNHcn\n0Te+w4AHShhfKa0xF2bWH/gD8DV3fzOBGEulzVy4+zvuvpm793H3zxP6Iwe7+xMJxVtM7XmP3A4M\nAch/4PUDFpY0ytJoTy5eIkxlj5ntAKxbxj0SY/VHwH8FjgHIHxG97fnb/7al5MNB7r7SzE4GphGK\n0JXuPs/MxgOPuvtdwJXAtWb2AvAm4R++7LQzFxcBnwFuzg+JvOTuhyQXdXG0Mxef+BPKdDioPblw\n96lmtp+ZzQVWAKeV49FyO18XpwF/NLMfEZrEx65+i9llZtcDtcDnzOxlYBywDuDufrm7321mB5jZ\nP4H3gG+3a7v504lERKQC6YphEZEKpiIgIlLBVARERCqYioCISAVTERARqWAqAiIiFUxFQESkgqkI\niIhUsP8H95VsZUlH2V4AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"dt = dx # experiment with this\n",
"\n",
"#current_f = f_advection\n",
"#current_f = f_heat\n",
"current_f = f_wave\n",
"\n",
"for i in range(5): # takes this many time steps at a time\n",
" current_u = rk4_step(current_u, current_t, dt, current_f)\n",
" current_t += dt\n",
"\n",
"pt.ylim([-0.25, 1.25])\n",
"pt.grid()\n",
"pt.plot(mesh, current_u[0])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}