# Fourth Order Runge-Kutta Method

1. Jul 20, 2009

### scarllatti

Hi
As a high school student I was working on an independent project and at last I found final differential equations that I want to integrate them numerically. I read a little bit and found out that Fourth Order Runge-Kutta Method is one of the good methods. By the way could you please help me to do that. I want to write a C++ program to use numerically integrated numbers to simulate something. So could please help me with that. I mean to integrate these four differential equations preferably in C++.
I give one of the four equations as a sample but all of them are in attachment so please refer to attachment to see them.

Û = -(K1 (sin z)^2 + K2 (cos z)^2 + (K1 - K2) * sin z * cos z u

Thanks very very much

#### Attached Files:

• ###### DE.jpg
File size:
23.2 KB
Views:
195
2. Jul 20, 2009

### rcgldr

http://en.wikipedia.org/wiki/Runge–Kutta_methods

Issues are covered here:

http://en.wikipedia.org/wiki/Stiff_equation

I described a simple interative trapezoidal predictor-corrector algorithm here:

also here:

wiki_predictor_corrector_proposed_example.htm

Since the trapezoidal algorithm is a linear approximation for each time step, the time steps still need to be very small for even the iterative algorithm to work. Wiki describes this as a predictor corrector method to implement the second method in this list:

The iterative process is needed for implict algorithms because the next step value is an input to it's own calculation:

yn+s = ... + f(yn+s, ...)

Wiki mentions this at the top of the multi-step article:

http://en.wikipedia.org/wiki/Multistep_method

Last edited: Jul 20, 2009
3. Jul 25, 2009

### MATLABdude

Re: Numerical Integration

I'd start with getting / borrowing a copy of "Numerical Recipes in C++"

However, the big question here is: do you know what an ODE is, and how you'd go about solving different classes of these by hand? Secondarily, do you have enough exposure to controls to understand what these equations signify?

If not, you may wish to reevaluate your independent project. If you do, hat's off to you!

4. Jul 25, 2009

### maverick_starstrider

Re: Numerical Integration

What you are trying to do is called "Finite element" methods/analysis. This is a numerical way of projecting a differential equation forward in time by treating its evolution as a number of small discrete steps. I have to say I am impressed that you are tackling such a thing in highschool. Ultimately, if you think as a mathematical derivative you think of $=\lim{\Delta \rightarrow 0}\frac{f(x+\Delta x}-f(x)}{x+\Delta x - x}$ however, instead of $\Delta x$ going to zero to get the deriviative, you let it equal a small number and get an approximation to the derivative (i.e. $f(x+\Delta x)=f(x)+\Delta x f'(x)$). So look up finite element methods and specifically what you want is the Verlet algorithm (it is the most accurate/numerically stable). and you will basically code the differential equation as is replacing the verlet equivalent with the derivatives as needed.

5. Jul 25, 2009

### D H

Staff Emeritus
No, he is not looking for finite element analysis. He is looking for techniques to numerically solve an ordinary differential equation initial value problem. Finite element techniques numerically solve a partial differential equation. There is a huge difference in complexity between ODEs and PDEs.

6. Jul 25, 2009

### maverick_starstrider

I guess I'm confused as to what the variables are in your DE's. Is there more than one?

7. Jul 25, 2009

### D H

Staff Emeritus
This is an ODE with four dependent variables: u, v, θ, and ω. There is but one independent variable: Time.

8. Jul 25, 2009

### scarllatti

Let me talk about that project maybe it can help you to really understand what I want to do. I want to mention this is not an school project or related to anything only I became interested in it so I wanted to do some works on it because I really love physics. I became interested in motion of falling papers and leaves so first I did many experiments to conclude some statistics data about fluttering and tumbling in their motion and how are they occur. After that I read some papers about visualizing it so I thought that it is a good idea to do it myself too. I mean I wanted to use the simplified equations that Tanabe and Kaneko used in their paper. So I want to use the equations from their paper. I know what ODEs are and I can solve basic ones too but I am not really mathematically comfortable with the final differential equations that I want to integrate them numerically as thier mentioned in their paper.
I attache the paper that I am specifically interested in to this reply. I will appreciate if someone can guide me through this process.
Thanks very much

#### Attached Files:

• ###### Behavior of a Falling Paper.pdf
File size:
733.1 KB
Views:
73
9. Jul 25, 2009

### D H

Staff Emeritus
Rather than starting with RK4, or any of the techniques Jeff mentioned in his post, you should start simple. The simplest technique for numerically solving ODEs is Euler's method. A couple of articles:
http://www.swarthmore.edu/NatSci/echeeve1/Ref/NumericInt/Euler1.html

Do a Google search "Euler's method" and you will find tons of online articles. The reason you need to understand Euler's method first is because all of the integration techniques mentioned use Euler's method as a starting point. Note well: Nobody uses Euler's method as an ending point is because it is, by itself, an incredibly lousy technique. Nonetheless, understanding Euler's method is critical for understanding more advanced techniques.

10. Jul 25, 2009

### rcgldr

The trapezoidal predictor corrector alogrithm improves Euler in two ways. First is using the average of the slopes at the current and next step. Second, iteration using the estimated next step value as an input provides feedback and will (normally) converge to a specific value, similar to Newtons root finding method (except the first guess is a given known good value):

initialize:
ynext = ycurrent

then iterate:
ynext = ycurrent + 1/2 (y'(ycurrent) + y'(ynext))

for some number of steps (6 to 16). ynext will converge after a reasonable number of steps if Δt is small enough.

The initialize and first iteration are essentially Euler

ynext = ycurrent + (y'(ycurrent))

but there's no need to code anything other than the trapezoidal (slope averaging) iteration step, after setting ynext = ycurrent.

4th order Runge-Kutta includes two weighted mid-point estimates, but Δt can just be halved to accomplish a similar thing. I doubt that complex algorithms designed to minimize computation time are actually needed much with the speeds of computers, including PC's these days.

I made an example Excel spreadsheet that shows this method used to solve y'=e-15y, with initial value y=1 at t=0, so y=e-15t. This is an example of a stiff equation that some numerical integration methods have issues with.

http://jeffareid.net/misc/numerical_integration.zip

A convergence test can be optionally added to the iteration loop from above:

iteration with convergence test:
yprev = ynext
ynext = ycurrent + 1/2 (y'(ycurrent) + y'(ynext))
if( | ynext - yprev | < tolerance ) early exit

Last edited: Jul 26, 2009
11. Aug 18, 2009

### CFDFEAGURU

I think one question that needs to be asked is "How accurate do you want to be?" MatlabDude set you onto the right website, www.nr.com (Numerical Recipes). Also, if you need to solve systems which most of the time you will since all higher order ODE's can be represented by a system of single order ODE's the implementation can be difficult and the NR guys already got it all figured out. (For the most part) Also, if the system is considered to be "Stiff" you will need much more advanced algoritms than RK4. For instance, Simi-implicit Euler. Also, you need to have the ability to automatically adapt the time stepping. This too is done by the NR guys.

Check out that resource, I live by it.

Thanks
Matt