{ "cells": [ { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Two versions of the leapfrog method is used to solve a system consisting of\n", "# a spring attached to a\"wall\" moving at velocity v. A mass m is attached to \n", "# the other end of the spring. The mass experiences friction as it moves (it\n", "# sits on a surface) and there is a rest friction as well that has to be overcome\n", "# in order for the mass to move at all.\n", "\n", "# The function accel() adds up all forces and returns the acceleration\n", "# There is a small inperfection in how the frictional force ffri is treated:\n", "# - initially when the mass is at rest, sign(v)=0 and there will be no friction\n", "# - at the fillowing step, the mass has velocity > 0, and ffri < 0\n", "# - the small v > 0 of the first step will the n rapidly decrease\n", "# - in some cases, a small negative velocity may result from the discretization\n", "# - but then the frictional force will be > 0 and again act to impede the motion\n", "# The result of the imperfection is that the mass initially, when it should\n", "# really stay at rest, has some very low, in some case oscillating, velocity\n", "# The behavior is \"self correcting\", diminishing when dt is decreased\n", "# - because of the self-correction there is no need to treat v=0 as a special case\n", "\n", "# The spring has native (relaxed) length lspr and can be elongated without bounds\n", "# and compressed to length 0. The final compression force is divergent to avoid the\n", "# spring length becoming negative (in principle a max length should also be imposed)\n", "# The friction force if fr0 at rest, decays rapidly to fr1 as v increases\n", "\n", " function accel(x,v,t,mass,kspr,lspr,xspr,kwall,fr0,fr1)\n", " fspr=(xspr-x-lspr)*kspr-kwall/abs(xspr-x)^6 # force from spring\n", " ffri=fr1+(fr0-fr1) # *exp(-abs(v)*10) # friction force magnitude\n", " ffri=-ffri*sign(v) # sign of friction depends on velocity\n", " return (fspr+ffri)/mass\n", " end\n", "\n", "# The function integrate1() carries out nt spes of the standard leapfrog algorithm\n", "# Results are written to the file \"x.dat\" every wt steps and is also pushed to vectors\n", "# The step error is not very good with the leapfrog method in the presence of friction\n", "# - good results can still be obtained quickly with reasonable dt values\n", "\n", " function integrate1(dt,wt,nt,mass,kspr,lspr,vspr,kwall,fr0,fr1)\n", " v=0.\n", " x=0.\n", " fi=open(\"x.dat\",\"w\")\n", " vect=Vector{Float64}()\n", " vecx=Vector{Float64}()\n", " vecv=Vector{Float64}()\n", " vecl=Vector{Float64}()\n", " for i=0:nt\n", " t=i*dt \n", " xspr=lspr+t*vspr \n", " if mod(i,wt)==0\n", " println(fi,t,\" \",x,\" \",xspr-x,\" \",v) \n", " push!(vect,t)\n", " push!(vecx,x)\n", " push!(vecv,v)\n", " push!(vecl,xspr-x)\n", " end\n", " v=v+dt*accel(x,v,t,mass,kspr,lspr,xspr,kwall,fr0,fr1)\n", " x=x+dt*v \n", " end\n", " close(fi)\n", " return vect,vecx,vecv,vecl\n", " end\n", "\n", "# The function integrate2() carries out nt steps of the leapfrog algorithm with\n", "# correction steps to better treat the damping. Otherwise same as integrate1()\n", "\n", " function integrate2(dt,wt,nt,mass,kspr,lspr,vspr,kwall,fr0,fr1)\n", " v=0.\n", " x=0.\n", " x1=0.\n", " fi=open(\"x.dat\",\"w\")\n", " vect=Vector{Float64}()\n", " vecx=Vector{Float64}()\n", " vecv=Vector{Float64}()\n", " vecl=Vector{Float64}()\n", " for i=0:nt\n", " t=i*dt \n", " xspr=lspr+t*vspr \n", " if mod(i,wt)==0\n", " println(fi,t,\" \",x,\" \",xspr-x,\" \",v) \n", " push!(vect,t)\n", " push!(vecx,x)\n", " push!(vecv,v)\n", " push!(vecl,xspr-x)\n", " end\n", " u=v+0.5*dt*accel(x,v,t,mass,kspr,lspr,xspr,kwall,fr0,fr1)\n", " y=x+dt*u # intermediate position\n", " u=0.5*(y-x1)/dt # intermediate velocity, x1 is x_{n-1}\n", " v=v+dt*accel(x,u,t,mass,kspr,lspr,xspr,kwall,fr0,fr1)\n", " x1=x # saving what will be x_{n-1} in the next step\n", " x=x+dt*v \n", " end\n", " close(fi)\n", " return vect,vecx,vecv,vecl\n", " end\n", "\n", "# Setting parameters related to the leapfrog loop\n", " dt=0.01 # time step\n", " wt=10 # write to file every wt step\n", " nt=5000 # number of time steps\n", "\n", "# Setting parameters related to the system \n", "# - mass connected to spring, the end of which is connected to a moving wall\n", "\n", " mass=1. # block mass\n", " kspr=1 # spring constant\n", " lspr=1. # native spring length\n", " vspr=0.5 # constant velocity of spring end (wall); a kind of driving\n", " kwall=0.0001 # constant in repulsive spring compression force\n", " fr0=4. # friction force at rest\n", " fr1=0.2 # friction force while in motion\n", "\n", "# Doing the integration for time step dt and then the same for dt/10,\n", "# plotting t,x,v, or l (the latter being the spring length) after each time.\n", "\n", " t,x,v,l=integrate1(dt,wt,nt,mass,kspr,lspr,vspr,kwall,fr0,fr1)\n", " \n", " using Plots\n", "\n", " plot(t,x)\n", "\n", " t,x,v,l=integrate2(dt,wt,nt,mass,kspr,lspr,vspr,kwall,fr0,fr1)\n", "\n", " plot!(t,x)\n", "\n", " dt=dt/10\n", " wt=wt*10\n", " nt=nt*10\n", "\n", " t,x,v,l=integrate1(dt,wt,nt,mass,kspr,lspr,vspr,kwall,fr0,fr1)\n", "\n", " plot!(t,x)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "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": 4 }