help with verlet integration, extended RK2, extended RK4
Hello everybody, I'm a beginner to c++ and I havent done any complex numerical calculations yet. I need some help with verlet integration.
I need to solve the following equation using verlet integration, extended RK2, and extended RK4:
y" = -g/L sin(y)
I have no idea where to start from. I would appreciate any help, thank you very much.
I need to solve the following equation using verlet integration, extended RK2, and extended RK4:
y" = -g/L sin(y)
I have no idea where to start from. I would appreciate any help, thank you very much.
Verlet integration
==================
Amounts to a standard approximation for the second derivative. Your equation is of the form
and the LHS is approximated to give
where the yi are successive values of your function. The subscripts i-1, i, i+1 etc just denote the time level.
Rearrange this as
Thus, you can get yi+1 from yi and yi-1 and hence march forward in time.
You probably don't need to store these values in an array, just have computational variables ynext, ycurrent and yprevious and output (to file and/or screen) successive values of ynext as you calculate and update these.
You will need two initial values of y; these depend on your boundary conditions (and what your professor has told you to do).
Note:-
(a) In your problem f = -(g/L) sin y
However, if you keep f as a separate function then you will easily be able to modify your code to solve other differential equations of similar form.
(b) I have used t (for time) as the independent variable, since your equation describes the angular motion of a pendulum with time. If you want to use x (or anything else) instead then that is fine.
(c) If y remains small then sin(y) is well approximated by y. Your equation then becomes simple harmonic motion (SHM), which can be solved analytically; you can thus check your answer. (In fact, if you let the function f return y instead of sin(y) that is precisely what you will get, so it would be a useful code-development strategy.)
Runge-Kutta (Your RK2 and RK4)
===========
I have no idea what your professor intends here, but I assume that "extended" refers to the fact that this is a second-order equation. Personally, I would rearrange it to a first-order equation as follows.
Rewrite as (sorry about the formatting, but I don't know any other way to get a column vector):
or
or, in vector form:
Lo and behold - a first-order differential equation! (but with vector dependent variable Y).
Now, your engineering maths textbook (or google) will give you the second-order (RK2) or fourth-order (RK4) Runge-Kutta method. If you just look for "Runge-Kutta" you will probably be given RK4: it's the commonest method in use in science and engineering. Matlab also has an RK4 solver to check your answers against.
Since Y is actually a two-component entity (first component is y, second is dy/dt) you can either:
(i) deal with each component separately;
(ii) put it in a 2-component array (actually, I would use a valarray in C++).
I would use a valarray since you can then write mathematical vector expressions like
A = B + C
which it will do componentwise automatically, just as you would write on paper. Standard Fortran arrays have done this by default since Fortran90. The nearest equivalent in C++ is a valarray.
This isn't a maths web forum: to go any further you need to post some of your code.
==================
Amounts to a standard approximation for the second derivative. Your equation is of the form
d2y/dt2 = f |
and the LHS is approximated to give
(yi+1- 2yi+yi-1)/dt2 = f(yi) |
where the yi are successive values of your function. The subscripts i-1, i, i+1 etc just denote the time level.
Rearrange this as
yi+1 = 2yi - yi-1 + dt2 f(yi) |
Thus, you can get yi+1 from yi and yi-1 and hence march forward in time.
You probably don't need to store these values in an array, just have computational variables ynext, ycurrent and yprevious and output (to file and/or screen) successive values of ynext as you calculate and update these.
You will need two initial values of y; these depend on your boundary conditions (and what your professor has told you to do).
Note:-
(a) In your problem f = -(g/L) sin y
However, if you keep f as a separate function then you will easily be able to modify your code to solve other differential equations of similar form.
(b) I have used t (for time) as the independent variable, since your equation describes the angular motion of a pendulum with time. If you want to use x (or anything else) instead then that is fine.
(c) If y remains small then sin(y) is well approximated by y. Your equation then becomes simple harmonic motion (SHM), which can be solved analytically; you can thus check your answer. (In fact, if you let the function f return y instead of sin(y) that is precisely what you will get, so it would be a useful code-development strategy.)
Runge-Kutta (Your RK2 and RK4)
===========
I have no idea what your professor intends here, but I assume that "extended" refers to the fact that this is a second-order equation. Personally, I would rearrange it to a first-order equation as follows.
d2y --- = f dt2 |
Rewrite as (sorry about the formatting, but I don't know any other way to get a column vector):
d ( y ) ( dy/dt ) ( dy/dt ) - ( ) = ( ) = ( ) dt( dy/dt) ( d2y/dt2) ( f ) |
or
d ( y0 ) ( y1 ) - ( ) = ( ) dt( y1 ) ( f ) |
or, in vector form:
dY -- = F dt |
Lo and behold - a first-order differential equation! (but with vector dependent variable Y).
Now, your engineering maths textbook (or google) will give you the second-order (RK2) or fourth-order (RK4) Runge-Kutta method. If you just look for "Runge-Kutta" you will probably be given RK4: it's the commonest method in use in science and engineering. Matlab also has an RK4 solver to check your answers against.
Since Y is actually a two-component entity (first component is y, second is dy/dt) you can either:
(i) deal with each component separately;
(ii) put it in a 2-component array (actually, I would use a valarray in C++).
I would use a valarray since you can then write mathematical vector expressions like
A = B + C
which it will do componentwise automatically, just as you would write on paper. Standard Fortran arrays have done this by default since Fortran90. The nearest equivalent in C++ is a valarray.
This isn't a maths web forum: to go any further you need to post some of your code.
Last edited on
Topic archived. No new replies allowed.