Transcript for:
Solving Differential Equations in Python

today we're solving differential equations in python this isn't necessarily an obvious or easy thing to do and there's a lot of subtle steps you need to take to prepare your problem to feed it into python properly and that's what we're looking at today we'll be looking at first order differential equations coupled first order differential equations how to prepare second order and higher order differential equations to pass into the psi pi functions that you use and even a few precautions you need to take as well when you're dealing with complicated or chaotic differential equations as always if you enjoy like and subscribe join the discord server there's a link in the description below to that and enjoy so most of the ode solving in python is done with scipy and you'll see that i have two important functions here for solving differential equations so other than that we're just going to be importing numpy and matplotlib so we'll start with first order odes because i think these establish a basis of what solving everything else is going to be like so i've chosen an ode here this is air friction when you're falling i've covered it in other videos be sure to check those out and so you have a dv dt minus alpha v squared this is just some constant and beta as a constant is equal to zero this is a differential equation with an initial condition the velocity at time zero is equal to zero so you're just standing still now if you want to solve an ode in python you want to make sure that you package everything in the correct form you can use those functions and that's one thing but you also want to package things properly as you get more and more into coding you realize that everything is about packaging stuff properly so we want to wrote write our ode in the form dv dt is equal to some function of t and v of course the differential equation can depend on time itself it doesn't in this case and it can also depend on that quantity that you're actually taking the derivative of so thinking about it in terms of english derivative of v equals something that depends on v and time that's what a differential equation really is in its essence so this equation is really simple we just move this over to the right hand side we get dv dt is equal to alpha v squared minus beta this that you see here is that function f of t v so we've written it in this form so we're starting to put our package together that we can feed into sci pi and like i said there's no dependence on time in this example so what we need to do is we need to write this in a python form and for that we need functions and so our function here we define a function dv dt right that's what you see here it's a function that we call dv dt and it depends on t and v and i always write it in this form because it helps to sort of keep things organized like this so you know before we deal with any ode solvers or anything we're just writing a function this is what the derivative is equal to it depends on time and v time comes first and we return okay well here alpha is equal to 3 and beta is equal to 5. so this is our differential equation it's written like this in a python functional form that we can eventually feed to a solver and i have the initial condition here of v of time equals zero is equal to zero so now that it's packaged in this form we can actually feed it into some of the scipy modules that will iterate through different algorithms these are very complex and they're sort of beyond the scope of this video but there's lots of different algorithms they use and you know those are dealt with by the solvers and i've sort of treated them like a black box in this video and so there's two main functions that scipy offers there's odent this is pretty classic and it uses a particular solver called um i think this is l soda it might be an iron l i don't remember and it's from a fortran library ode pack so it's a particular solver when you use odent alternatively there's solve ivp this is far more customizable there's a bunch of different solvers you can use and it's actually a parameter to the solver but when you use these two functions you sort of pass them in the same things there's a little bit of differences and they spit out something that's a little bit different so we have to get used to each so here's my differential equation here and so here i say okay i want to solve for a number of times i'm going to solve for a hundred times between zero and one and i wanna plot the solution between zero and one second 100 times that's what linspace does of course if i just copy and paste this i believe i need to run all the code as well so here's my time these are the times at which i want to get the solution you can see here i want to find what v is and so i'm going to solve it using both ode solvers i have ode end here it takes in this differential equation remember that this is the function that we defined above dv dt it takes in our initial condition it takes it as a just a scalar quantity v naught so why not y naught is the name of the parameter in this function and our initial condition is v naught so it calls it y naught the argument in the function t is the times we want to solve it over so t equals t that's just the times we define here and t first is equal to true this is because odent and solve ivp are slightly different you'll note that in this function i have t comes first otherwise odent actually expects v to come first so i like to write it the same way because it'll work for solve ivp or odent so i'm going to specify that t first is equal to true this will give us one solution and you'll see it returns an interesting object uh solve ivp same thing we give it the function dv dt then we have to give it a t span so we it needs to know how long am i solving this over and it's solving between zero and whatever the maximum value of this array is which is one and y naught is equal to v naught you'll note that i give it in an array form this is seems weird when you're solving one dimensional odes or just a system of one ode it makes more sense when you're solving systems of odes and that's really what python is for is for systems of coupled differential equations and then the initial condition of course is just an array of many many different things and then the times to evaluate on t val is t that's this array here so we can solve our ode and we can look at the two different solutions and they return slightly different objects so solution m1 okay it looks weird it's a 2d array and we have one thing each here this is because if there were multiple quantities you were solving for which we'll get to in a bit if you have you know multiple you know maybe you're solving for y1 y2 y3 you know coupled things together it would return like y1 then y2 then y3 and this corresponds to the first time the second time the third time so because we're only solving for one quantity here which is velocity uh it only returns one thing but if there were multiple things you know this would be maybe two values and this would be two values as well and so if i want to extract a quantity quick i can go transpose this array right and i take the first element and that gives me my velocity now obviously there's nothing else but if there were multiple things in this um differential equation or if it was coupled or we had multiple differential equations together then i can index to get the other things as well so this is how i would get my velocity as a function of time like this uh we could also look at solution m2 this is a little bit different how it returns from solve ivp it actually returns a more complicated object where it says okay the solver successfully found the solution it gives the number of iterations and a bunch of stuff it says okay it's successful and it gives this t which are the times that it solves for and it also gives t events which is something else and then it gives y and y is the solution to the ode so if i want the solution i have to go dot y and it will give y and like before um you'll see that it's a 2d array and if there were multiple quantities it would list them sort of um you know maybe this would be y1 then there'd be another thing after y2 y3 etc etc so they return slightly different objects it's use getting it's good to get used to both of them and typically i'd probably use solve ivp more often just because of the versatility of the solver so here i'm extracting both the solutions this is the v solution so the velocity using method one method one is ode int and i just extract it like i did before and solution m2 i extract the velocity here using exactly what i did before as well note that this is transpose i transpose an array i take the first element and this is of course getting the solution in that weird solve ivp object so i get both velocities and i can plot them and you'll see that it indeed solves the differential equation the same way using both solvers so this is an introduction to both solvers here so what python is really useful for is coupled first order differential equations and these are really the things that will show up all the time as you get further and further into research going into topics et cetera and so this is an example of a couple first order differential equation we have y1 and y2 these are two quantities that we want to solve for and y1 prime is the derivative of y1 and y2 prime is the derivative of y2 but the derivative of y1 also depends on y2 and the derivative of y2 also depends on y1 furthermore there's this also dependence on x so x is kind of like t in this example and i purposefully chose to use x instead of t because often you'll be dealing with differential equations with variables that don't match right like t here the arguments say t here but you have to think of t as an x so it's used getting used to this is like the independent variable of the differential equation the thing that you change and how everything else changes so when you're solving coupled first order odes there's a strategy that i always use and i would always recommend to do this because otherwise it gets sort of confusing so there are two quantities that we want we want y1 and y2 right and there are two initial conditions for y1 and y2 so i define a vector capital s and you'll see in all the videos that i do i always use this notation of capital s s is just a vector of everything that you want to solve for it's very easy so this is a vector that contains y 1 and y 2. before we pass our differential equations to an ode solver we have to define a function dsdx and of course that's just the derivative of these two things the function dsdx takes an s so dsdx which is a vector depends on s which is a vector and x right so looking at it like this so here's our vector s is equal to y1 y2 ds dx is a function it's a vector function of x and s so you can think of that as sort of analogous to this here dv dt is a function of t and v now it's in a vector form because we have um multiple variables ds dx is equal to f of x and s so this is like time and s is the uh the thing here and you'll note that this can be written like this right so s is a vector with y1 and y2 so this is really a function of that depends on x y1 and y2 and of course ds dx is just the derivative of y1 and y2 so y1 prime y2 prime which is equal to everything on the right hand side here so you'll note that what it i have above here everything is contained down here like this in ds dx so what i do in my python function is i write a function ds dx and it takes an x and it takes an s exactly like i have here now s is equal to y1 y2 so i extract those components from s and it returns dsdx so it returns whatever is here the right hand side of these differential equations so all the information about the differential equations is contained here so this is y one prime here and y two prime so my dsdx is simple like this i have the initial conditions y one and y two uh that i specify here of course those can be numbers or whatever that's fine and uh i put it in initial condition s naught so that's like the initial condition so you can think of s as like a vector that contains all the information i want for the system so i specify this and we can solve this ode so i want to solve for a bunch of values between x's 0 and 1 and i'm going to use ode interior so i pass my dsdx function exactly that's sort of the same like above my initial condition why not this is the argument name in the function and i pass an s naught which is a vector here t which are the times i want to solve for remember it's not times here this is equal to x which are the x's that i want to solve over i'm also using t first is equal to true so i can solve this and i can look at my solution and you'll see what i mean so these two values correspond to the first time these two values correspond to the second time the third time the fourth time are well x here because x is an array so this is the first time the second time the third time fourth time everything is discrete because you're on a computer of course so this is the first time the second time the third time the fourth time the fifth time this is y one y two at the first time all right first x uh y one y two at the second x y one y two and the third x so it's used getting you got to think of it sort of like in this way if i want to extract like a column of course what i would do is i would transpose this and take the 0th component that will give me y1 the next component would give me y2 so very simple right that's what i do here i get y1 and y2 from my solution like this and i can plot them and so what i've done is i've solved my differential equation and this is what the solution for y1 and y2 look like so you might be thinking everything we've done so far is for first order differential equations what if i want to do higher order differential equations like second order third order fourth order or anything like that well you'll see that any higher order differential equation can be converted into a system of first order differential equations so python doesn't have functions to directly solve second order odes you need to take your second order differential equation convert it into two first order differential equations if you had a third order differential equation that would be three first order or fourth order differential equation would be three first order differential equations and it keeps going sort of up like that and i'll show you what you mean so any second order ode can be verted can can be converted into two first order odes and i'll show you how to do that here so consider the differential equation x dot dot is equal to minus x dot squared plus sine x uh this is a function of time for example so derivative of a dot means with respect to time so i want to convert this into two first order odes so i can solve them like i did in the example above with two first order equations uh so what i do is i take x this is what i want to solve for and i define x dot is equal to v right so v is a new variable here this is a trick and you get used to this trick over time the first time you see it it might be a little confusing so v is our new variable and x dot equals v is a differential equation right um it's weird i've defined a new thing but i'm saying the derivative of x dot it depends on v so okay that it's just a definition but it is a differential equation and the other differential equation is v dot is equal to minus x dot squared plus sine x and x dot of course is just v so this is minus v squared plus sine x so we get two equations here just by introducing this definition this becomes one of our differential equations and our second differential equation which is with v v dot is equal to minus v squared plus sine x so we're solving for both v and x and so there are two coupled first order equations so like before i define this ds dx where of course s contains everything we want to solve for so s is a vector of x and v and it returns the derivatives here i specify an initial condition for x and v and s naught the initial condition vector and i can solve it like i did before so everything is the same and i get my solution here so the important thing about this part is you have a higher order differential equation you convert it into two first order differential equations using this trick here then you can package everything properly to put it into a python solver to really show off this here and especially for high order differential equations i'll give an example that sort of summarizes everything so far so here i have two coupled third order equations with x one and x two so i have a triple dot here it's like a third derivative and of course these functions are a little bit weird so x two dot dot dot and x one dot dot dot and i need to convert this into many first order equations so for that i need to do a bunch of definitions i'll define four new variables here v1 is x dot v2 is x2 dot a1 is x1 dot dot and a2 is x2 dot dot you'll note that this here which is dot dot dot i need to define new variables up to the dot before that so if i have x two like the third derivative i'll need to define a new variable that's the second derivative if this was like five dots i would need to define all the way up to a new variable that's like four dots so it's one less than the number of dots here that's how i'm defining the variable and then everything leading up to that as well so i define four new things here and these are four differential equations like before these are simple i'm just introducing new variables but they signify first order differential equations finally um a1 and a2 which are the dot before the one that comes here if i take a one dot that's x one dot dot dot right that's why i define one less than the one here because then i can take the derivative of this one and i get the quantity here and then it's just everything on this side so uh a1 dot is equal to x1 dot dot dot which is this here and a2 dot is x2 dot dot dot which is this here and then i just plug in everything from up here so a2 dot is equal to this it's just a1 cubed of course a1 is x1 dot dot so a1 cubed plus v2 that's x2 dot plus x1 plus sine t and these make up our final two differential equations so we have six coupled first order differential equations here we've turned two coupled third order equations here into six coupled first order differential equations i encourage you to think about this and go through it really slowly because this is essential for all differential equations solving in python so these are our six equations and now i define s so s is everything we want to solve for we want to solve for x1 and x2 but we also have these new variables that we defined v1 v2 and a1 and a2 so s is x1 v1 a1 x2 v2 a2 these are all six things that we want to solve for then dsdt is just of course the derivative for all these things so x1 dot is equal to v1 that's by definition v1 dot is a1 that's by definition a1 dot is minus v2 2v2 squared plus x2 that's just this equation here essentially same thing with x2 v2 and a2 on this equation here so my ds dt like before d it's a function that takes in s and t so i define ds dt it takes an s and t s is just all these quantities here x1 v1 a1 x2 v2 a2 and it returns the derivatives listed exactly like i have them in this column vector here i specify all the initial conditions of course you can change these and my initial condition vector here i run this and like before i just solve exactly the same i pass the dsdt my initial condition the times i want to solve for and making sure t first is true because t comes first in this function so i can solve it and i get a little solution like this and this isn't any like significant differential equation or anything but you see more importantly the process so if you run into some math that you're doing you get a complicated differential equation you know exactly how to set it up in python to solve it i'll just give a final note here and that final note is to be careful when you're solving odes not all solvers work for all odes if you use these solvers like a black box sometimes you'll get the wrong solution so if you're doing research and you get an ode make sure you know everything about that differential equation so in one of my videos and there's a link here and i'll put a link in the description to three body planetary motion you get a system of i don't know how many of these are your like 12 differential equations 12 coupled first order equations so these are system of odes but you can't just use an ode solver out of the box it will fail and it will diverge and it'll give some crazy solution that's not the right solution it requires the dope 853 solver that's you can use that in scipy and it also requires low values for arthole and atol these are parameters for the solve ivp method and there's a paper here about this system of differential equations and so you always need to review the literature about your specific ode especially if it's a complicated one before trying to solve in psi pi because otherwise you might get a solution that's just chaotic and so um you need to play around with artol and atl like i said so here's dsdt and i'm going to put in all these equations here so s contains all these variables and i define some things this is just for the denominator here because it gets a little um annoying to specify this every time and uh you know here's the quantities this is the sdt that it's returning and i define my initial conditions so there's 12 different things there's 12 initial conditions and i want to solve this so here's my t and here's solving it properly right and so if i solve it properly and i look at y1 for example which is the y position of the first mass i get something like this and this is the gravitational body sort of moving around each other but what happens if i don't set these parameters here that are required well okay i can try solving it and you'll note i get a completely different solution and this is small perturbation since it's a chaotic system if the solver does one thing slightly incorrectly it can lead to a sort of chaotic and solution that doesn't really make any sense and so we need to make sure we specify those are toll and atoll same thing if i don't use this method dope 853 if i just use a regular solver it'll run into another sort of chaotic solution so if you look at the pattern here this is one bad solution i plot it again and there's another bad solution so make sure you have the right solver and review the literature for that and also play with atoll and artol and these are like the tolerance for the ode solver at each step when i get new values for my solution how careful do i need to be of course the smaller those parameters are the longer it takes to solve the ode but setting them small gives a better or more precise result for the ode and it's essential for chaotic differential equations as well anyways i hope you enjoyed this video there's plenty more content if you like this video i have tons of videos in physics where i solve different physics problems that all involve complicated coupled ordinary differential equations so look at my channel any of the lagrangian mechanics double pendulum problems the three-body gravitational problem anything like that is all solving systems of coupled differential equations be sure to like subscribe join the discord server check out all the other videos and i'll see you next time