VIEWS: 4 PAGES: 4 POSTED ON: 5/25/2012
Solving a 1st Order Ordinary Differential Equation in Excel The Falling Ball Example Let us take the falling ball example and develop a numerical solution to the differential equation that describes its motion. We began with the differential equation: my mg FD where the positive y direction is down, m is the mass of the ball, and FD is the drag force exerted by the air. This is a second order ordinary differential equation, and although we will learn how to solve for the position of the ball as a function of time, y(t), in the near future, today we will just solve for the velocity of the ball as a function of time v(t). Stating our differential equation in terms of velocity, and dividing by the mass, we have: v g (1 FD / m) Actually, as we saw from our derivation in class, using the following definition for the terminal velocity, 1/ 2 mg vt 2 # R We can restate our differential equation as: v 2 v 9.8 1 vt Here we have written g as 9.8 m/s. To solve our differential equation, we must enter this equation for the derivative of the velocity into our Excel spreadsheet. To begin, start up Excel. Then go to the File menu, and then click on New… Select Workbook. A blank spreadsheet should appear in front of you. Let us now define a funtion representing the differential equation we have been discussing. Go to the Tools menu, select Macro and then Visual Basic Editor. A new window will open, placing you in the Visual Basic environment. We will be using this enviroment in a fairly simple way to create user defined functions for Excel. Select Module, then Insert. We are now ready to enter our function into the window that is probably titled “Book 1 – Module 1 (Code)”. To enter the function, type: Function Fun(t, v, vt) Fun = 9.8 * (1 - (v / vt) ^ 2) End Function As one can see, the function is entered as a set of basic commands. The first line is the header of the function. Functions must always start with the keyword “Function” followed by the name of the function as you will refer to it in the Excel spreadsheet. Following the name of the function, “Fun” (in this case) is a list of parameters that the Excel spreadsheet will pass to the function. We will pass the time, t, the current velocity, v, and the terminal velocity, vt. On the next line we assign the value of the function according to our differential equation. This particular function represents the derivative of the velocity for our falling ball. On the following line, we signal the end of the function with the keywords “End Function”. Go to the File menu of the Visual Basic Editor and select Close and Return to Microsoft Excel. Now we are ready to use the 4th order Runge Kutta method to solve our differential equation. In fact, let us compare the 4th Order Runge-Kutta method with both the Euler method and the analytical solution. You will need to define some constants. First create a series of labels: grav 9.8 m/s roair 1.21 kg/m^3 roball 1200 kg/m^3 radius 0.04 m To define the grav constant, select the cell containing the number 9.8. Then select Insert, Name, and Define. Excel is smart enough to guess that you want to call the value 9.8 “grav”, so just click okay. Do the same for the other constants. From these, you can calculate and then define the derived constants: mass 0.321699 Kg vt 22.76724 m/s Down about 10 rows I enter the headings for my Excel Spreadsheet: grav 9.8 m/s roair 1.21 kg/m^3 roball 1200 kg/m^3 radius 0.04 M mass 0.321699 Kg vt 22.76724 m/s delta 0.5 S Analytical Euler's Method 4th Order Runge-Kutta Method for 1st Order ODE time velocity Accel velocity accel velocity f k1 k2 k3 k4 For the time column, enter 0 for the first time. Let’s say the 0 is entered in cell A23. Go down one cell further to A24, and type: = A23 + delta. You can then fill down as many cells as you like. Fill down so that times 0 to 20 seconds are covered. The anaytical solution for velocity is: =vt*(EXP(2*grav*A23/vt)-1)/(EXP(2*grav*A23/vt)+1) Where A23 refers to the time. The analytical solution for the acceleration is: =grav*(1-(B23/vt)^2) This is just our original diffential equation for acceleration, with our Solution for the velocity, B23, plugged in. Our first timestep for Euler’s method velocity is the initial velocity. Use 0. Our first acceleration for the Euler’s method is =grav*(1-(D23/vt)^2) where D23 is the Euler’s method velocity for that time step. The second timestep for Euler’s method velocity is calculated: =D23+delta*E23 That is, the new velocity is the old velocity plus the timestep length times the accleration. As for the Runge-Kutta Columns, first enter 0 for the intial velocity. The column labeled f calculates the acceleration using the current velocity. =Fun(A23,F23,vt) We will discuss the remaining columns in class.