I'm trying to simulate a chain of coupled oscillators and to do so I need to use Runge-Kutta order 4 with an array (for all the different oscillators). My problem is that I'm really bad at declaring variables correctly and RK4 isn't doing what it should probably because I'm using a pointer for the calculation. (using one because it won't work without it idk why). if anyone knows how I can work my way around this problem would be awesome. Thank you!
My functions f1 and f2 are in a header file. Any function that doesn't depend on my index i work perfectly well. (Typically if f2 returns 1, I get as expected 1/2x^2 for the second integral). As soon as I put f2 that depends on i I always get something that diverges. Code below: (good luck)
HEADER FILE:
#include <iostream>
#include <cmath>
using namespace std;
double f1(double theta1) {
return theta1;
}
double f2(double t, double theta[], double alpha, double beta, int i) {
//return sin(theta[i])-beta*(theta[i+1]+theta[i-1]-2*theta[i]);
return sin(theta[i]);
}
MAIN FILE:
#include "func_solition.h"
#include <iostream>
using namespace std;
int main() {
// Paramètres intégrale
double dt = 0.01;
double t = 0;
double TF = 20;
// CI pendules
int n = 10; // ATTENTION ON A DEUX PENDULES POUR LES BC, PREVOIR +2 PENDULES
double theta[n];
double theta1[n];
for (int i = 0; i < n; i++) {
theta[i] = 0;
}
theta[1] = 3.141592/2;
for (int i = 0; i < n; i++) {
theta1[i] = 0;
}
theta1[1] = 0;
// Constantes équation du mouvement
double alpha = 0.1;
double beta = 0.01;
// Variables pour rk4
double tmp1=0.0, tmp2=0.0;
double k1, k2, k3, k4;
double l1, l2, l3, l4;
cout << "CI: theta[1] = " << theta[1] << ", theta1[1] = " << theta1[1] << endl;
// rk4
while (t <= TF) {
for (int i = 1; i < n-2; ++i) {
k1 = dt * f1(theta1[i]);
l1 = dt * f2(t, theta, alpha, beta, i);
k2 = dt * f1(theta1[i] + 0.5 * l1);
tmp2 = theta[i] + 0.5 * k1;
tmp1 = t + 0.5 * dt;
l2 = dt * f2(tmp1, &tmp2, alpha, beta, i);
k3 = dt * f1(theta1[i] + 0.5 * l2);
tmp2 = theta[i] + 0.5 * k2;
tmp1 = t + 0.5 * dt;
l3 = dt * f2(tmp1, &tmp2, alpha, beta, i);
k4 = dt * f1(theta1[i] + l3);
tmp2 = theta[i] + k3;
tmp1 = t + dt;
l4 = dt * f2(tmp1, &tmp2, alpha, beta, i);
theta[i] = theta[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
theta1[i] = theta1[i] + (l1 + 2 * l2 + 2 * l3 + l4) / 6;
}
> Boundary conditions
theta[n-1] = 0;
theta1[n-1] = 0;
theta[0] = 0;
theta1[0] = 0;
t = t + dt;
cout << t << ", " << theta[2] << ", " << theta1[2] << endl;
}
}
I'm trying to simulate a chain of coupled oscillators and to do so I need to use Runge-Kutta order 4 with an array (for all the different oscillators). My problem is that I'm really bad at declaring variables correctly and RK4 isn't doing what it should probably because I'm using a pointer for the calculation. (using one because it won't work without it idk why). if anyone knows how I can work my way around this problem would be awesome. Thank you!
My functions f1 and f2 are in a header file. Any function that doesn't depend on my index i work perfectly well. (Typically if f2 returns 1, I get as expected 1/2x^2 for the second integral). As soon as I put f2 that depends on i I always get something that diverges. Code below: (good luck)
HEADER FILE:
#include <iostream>
#include <cmath>
using namespace std;
double f1(double theta1) {
return theta1;
}
double f2(double t, double theta[], double alpha, double beta, int i) {
//return sin(theta[i])-beta*(theta[i+1]+theta[i-1]-2*theta[i]);
return sin(theta[i]);
}
MAIN FILE:
#include "func_solition.h"
#include <iostream>
using namespace std;
int main() {
// Paramètres intégrale
double dt = 0.01;
double t = 0;
double TF = 20;
// CI pendules
int n = 10; // ATTENTION ON A DEUX PENDULES POUR LES BC, PREVOIR +2 PENDULES
double theta[n];
double theta1[n];
for (int i = 0; i < n; i++) {
theta[i] = 0;
}
theta[1] = 3.141592/2;
for (int i = 0; i < n; i++) {
theta1[i] = 0;
}
theta1[1] = 0;
// Constantes équation du mouvement
double alpha = 0.1;
double beta = 0.01;
// Variables pour rk4
double tmp1=0.0, tmp2=0.0;
double k1, k2, k3, k4;
double l1, l2, l3, l4;
cout << "CI: theta[1] = " << theta[1] << ", theta1[1] = " << theta1[1] << endl;
// rk4
while (t <= TF) {
for (int i = 1; i < n-2; ++i) {
k1 = dt * f1(theta1[i]);
l1 = dt * f2(t, theta, alpha, beta, i);
k2 = dt * f1(theta1[i] + 0.5 * l1);
tmp2 = theta[i] + 0.5 * k1;
tmp1 = t + 0.5 * dt;
l2 = dt * f2(tmp1, &tmp2, alpha, beta, i);
k3 = dt * f1(theta1[i] + 0.5 * l2);
tmp2 = theta[i] + 0.5 * k2;
tmp1 = t + 0.5 * dt;
l3 = dt * f2(tmp1, &tmp2, alpha, beta, i);
k4 = dt * f1(theta1[i] + l3);
tmp2 = theta[i] + k3;
tmp1 = t + dt;
l4 = dt * f2(tmp1, &tmp2, alpha, beta, i);
theta[i] = theta[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
theta1[i] = theta1[i] + (l1 + 2 * l2 + 2 * l3 + l4) / 6;
}
> Boundary conditions
theta[n-1] = 0;
theta1[n-1] = 0;
theta[0] = 0;
theta1[0] = 0;
t = t + dt;
cout << t << ", " << theta[2] << ", " << theta1[2] << endl;
}
}
Share
Improve this question
edited Feb 17 at 15:42
3CxEZiVlQ
38.8k11 gold badges76 silver badges90 bronze badges
asked Feb 17 at 15:39
user29681619user29681619
211 bronze badge
New contributor
user29681619 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
12
|
Show 7 more comments
1 Answer
Reset to default 2You are having to loop through your theta (and theta derivative) arrays. However, if you use a mathematics-understanding array like valarray<double>
then you can treat the array as a single variable and use the simple RK4 formula.
Note that both angle (theta) and its derivative (theta1) are contained in the same array: angles first, then angular velocity. You need a single routine to return the derivative of the array. For the latter half of the array it will involve your angular-acceleration function.
In the code below I'm guessing a bit what you really intended as an angular acceleration function. Note: I have increased your value of beta quite a bit, or the output isn't very interesting!
To see its output, redirect the output to a file:
a.exe > out
and then plot with gnuplot:
p for [i = 2:9] 'out' u 1:i+1 w l title sprintf("theta[%d]",i)
Code:
#include <iostream>
#include <valarray>
#include <cmath>
using namespace std;
using vec = valarray<double>;
//======================================================================
vec derivative( double t, const vec &theta )
{
const double alpha = 0.1;
const double beta = 0.1;
int n2 = theta.size(); // angles, then angular velocity
int n = n2 / 2; // number of pendula
vec result( n2 ); // initialised to 0 by default
for ( int i = 1; i < n - 1; i++ ) // derivative of angle is a later part of the vector
{
result[i] = theta[i + n];
result[i+n] = -alpha * sin( theta[i] ) + beta * ( theta[i-1] - 2 * theta[i] + theta[i+1] );
}
return result;
}
//======================================================================
vec rk4step( double t, const vec &Y, vec f( double, const vec& ), double dt )
{
vec dY1, dY2, dY3, dY4;
dY1 = dt * f( t , Y );
dY2 = dt * f( t + 0.5 * dt, Y + 0.5 * dY1 );
dY3 = dt * f( t + 0.5 * dt, Y + 0.5 * dY2 );
dY4 = dt * f( t + dt, Y + dY3 );
return ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;
}
//======================================================================
int main()
{
// Timestepping parameters
double dt = 0.01;
double tmax = 20;
// Pendulum parameters
int n = 10; // number of pendula (first and last do not oscillate)
vec theta( 2 * n ); // n angles, followed by n angular velocities
theta[1] = 3.141592 / 2; // draw the first pendulum aside and then release ...
double t = 0;
while ( t <= tmax )
{
theta += rk4step( t, theta, derivative, dt );
t += dt;
cout << t;
for ( int i = 0; i < n; i++ ) cout << " " << theta[i];
cout << '\n';
}
}
Output:
f2
expects an array for 2nd parametertheta
, but you pass only one number pretending to be an array -&tmp2
. That cannot possibly end well. May I suggest getting a good C++ book to learn from before you delve into implementing complicated maths algorithms? – Yksisarvinen Commented Feb 17 at 15:42std::array
orstd::vector
? – Jesper Juhl Commented Feb 17 at 15:49