Hello,
Thanks for taking the time to reply! However, I had already solved the problem a few days ago on my own. The DOPRI 8(7)13 is the Dormand Prince 8th order with 13 stages method for solving <http://en.wikipedia.org/wiki/Ordinary_differential_equations> ordinary differential equations. It is a member of the Runge-Kutta family of ODE solvers. I have been able to reproduce it in Python with a multi-variable capability. On a particular test equation the error of the RK4 is about 1e-4 whereas the DOPRI 8(7)13 is about 1e-13. I have then used the multi-variable DOPRI to compute for the n-body problem using 1st order ODEs which requires you to have 4 equations for 1 planet. I have already made a 10-body problem through an n-body solver in python and using DOPRI as my integrator. It works fine. Thanks again for your time and consideration. Regards Vick -----Original Message----- From: Oscar Benjamin [mailto:oscar.j.benja...@gmail.com] Sent: Monday, 22 July, 2013 15:59 To: Vick Cc: Tutor@python.org Subject: Re: [Tutor] hi On 12 July 2013 11:08, Vick < <mailto:vick1...@orange.mu> vick1...@orange.mu> wrote: > Hi, Hi Vick, Sorry for the delayed response I didn't see this email until a while after you sent it. > I have written a code to perform a numerical solution to 1st order > Ordinary Differential Equations (ODEs). I have written codes for the > famous Runge Kutta 4th order and the DOPRI 8(7)13. However the codes > are for 1 variable only; viz. the "y" variable. It would be good to show us the code for this in your email (not in an attachment). Also the code should be complete: the code uses a function rungkdp8 which is not shown. Where did that come from? Until you know what you're doing you should test your integrator on a simpler set of equations. Also when posting to a mailing list you should simplify the problem as much as possible which would mean using simpler ODEs. In short read this: <http://sscce.org/> http://sscce.org/ The code attached is incorrect as it applies each stage of the Runge-Kutta method to all time-steps sequentially before moving on to the next stage. It should apply each stage for one time step and then use the output of that for the next time-step. In other words you need to change the order of the loops. Here is an example of how to implement an Euler stepper that should hopefully get you started on the right track. I've deliberately avoided using numpy in this example even though it makes most of this easier. I did this to focus on the basics of Python which you should learn first. # integrator.py nan = float('nan') # indices for the two variables in the system POS = 0 VEL = 1 VARS = ['pos', 'vel'] # system has one parameter OMEGA OMEGA = 10 # Frequency in Hz is OMEGA / 2pi # function for the derivative of the system def shm(t, x, dxdt): '''Equation of motion for a mass in a spring with freq OMEGA''' dxdt[POS] = x[VEL] # Python uses square brackets for indexing dxdt[VEL] = - OMEGA**2 * x[POS] # Initial conditions for the problem x0 = [nan] * 2 x0[POS] = 1 x0[VEL] = 0 # This is the function that you should replace with e.g. rk4 or dopri8 def euler(t1, x1, t2): '''Take 1 Euler step from x1 at t1 to x2 at t2''' # Create empty lists for answer x2 = [nan] * len(x1) dxdt = [nan] * len(x1) # Call derivative function to fill in dxdt shm(t1, x1, dxdt) # Compute x2 and return dt = t2 - t1 for n in range(len(x1)): x2[n] = x1[n] + dt * dxdt[n] return x2 def solve(ts, x0): '''Compute states corresponding to the times in ts x0 is the state at ts[0] (the initial condition). ''' # Create an empty list of lists for the solution Xt = [x0[:]] # The initial conditions # Loop through timesteps computing the next value for n in range(len(ts) - 1): Xt.append(euler(ts[n], Xt[n], ts[n+1])) return Xt def print_line(*items): print(', '.join(str(i) for i in items)) def print_solution(ts, Xt, vars): print_line('t', *vars) for t, x in zip(times, Xt): print_line(t, *x) # Numerical parameters for the solution DT = 0.001 t0 = 0 T = 1 times = [t0] while times[-1] < T: times.append(times[-1] + DT) # Solve and print Xt = solve(times, x0) print_solution(times, Xt, VARS) > I have written another code in Excel VBA for DOPRI 8 for a > multi-variable capability. However I have not been able to reproduce > it in Python. I'm having trouble in making arrays or lists, I don't > know which is supposed to work. > > I have attached my Excel VBA code for a multi-variable numerical > integrator in PDF format. This does work in Excel. I have also attached my python code. I'm not familiar with Excel VBA but I think that this code suffers from the same problem that the loop over stages takes place outside the loop over time-steps so I believe that it is incorrect. > Can anyone help me out please? I happen to know quite a lot about this subject so feel free to ask more questions. Oscar
_______________________________________________ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: http://mail.python.org/mailman/listinfo/tutor
