C++

Forum

 
How do I fix segmentation faults with vector iterations?

fusion809 (3)
Hi,

I would like to solve the ODE:

d2x/dt2=-9.8*cos(x)

x(0) = dx/dt(0) = 0

on t in [0,10], using the Runge-Kutta 4th order method (RK4). I would like the data the comes out to be stored in vectors called t, x and dx. Each having N+1 (or 1,001 in my code) elements and the step size being 10/1000 = 0.01. Here is my code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
 
#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <cmath>
#include <utility>

using namespace std;

double dx2(int t, int x, int dx)
{
	return (-9.8*cos(x));
}

std::pair<double, double> RK4(float t, float x, float dx, float h)
{
	double k1, k2, k3, k4, l1, l2, l3, l4, diff1, diff2;
	k1 = h*dx2(t,x,dx);
	l1 = h*k1;
	k2 = h*dx2(t+h/2,x+l1/2,dx+k1/2);
	l2 = h*k2;
	k3 = h*dx2(t+h/2,x+l2/2,dx+k2/2);
	l3 = h*k3;
	k4 = h*dx2(t+h,x+l3,dx+k3);
	l4 = h*k4;
    diff1 = (l1+2*l2+2*l3+l4)/float(6);
    diff2 = (k1+2*k2+2*k3+k4)/float(6);
    return {diff1, diff2};
}

int main()
{
    double t0, t1, x0, dx0, h, N;
    N     = 1000;
    t0    = 0;
    std::vector<int> t;
    t[0]  = t0;
    t1    = 10;
    x0    = 0;
    std::vector<int> x;
    x[0]  = x0;
    dx0   = 0;
    std::vector<int> dx;
    dx[0] = dx0;
    h     = (t1 - t0) / float(N);
    
    for(int i = 1; i<=N; i++) {
        auto diff = RK4(t[i-1],x[i-1],dx[i-1],h);
        x.push_back(  x[i-1] + diff.first  );
        dx.push_back(dx[i-1] + diff.second );
        t.push_back( t[i-1]  + h);
    }
    cout << x[N];
    return 0;
}


I keep getting segmentation faults whenever I execute this code. I have tried using this alternate code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
 
#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <cmath>
#include <utility>

using namespace std;

double dx2(int t, int x, int dx)
{
	return (-9.8*cos(x));
}

std::pair<double, double> RK4(float t, float x, float dx, float h)
{
	double k1, k2, k3, k4, l1, l2, l3, l4, diff1, diff2;
	k1 = h*dx2(t,x,dx);
	l1 = h*k1;
	k2 = h*dx2(t+h/2,x+l1/2,dx+k1/2);
	l2 = h*k2;
	k3 = h*dx2(t+h/2,x+l2/2,dx+k2/2);
	l3 = h*k3;
	k4 = h*dx2(t+h,x+l3,dx+k3);
	l4 = h*k4;
    diff1 = (l1+2*l2+2*l3+l4)/float(6);
    diff2 = (k1+2*k2+2*k3+k4)/float(6);
    return {diff1, diff2};
}

int main()
{
    double t0, t1, x0, dx0, h, N;
    N     = 1000;
    t0    = 0;
    int myints[] = {t0};
	std::vector<int> t(myints, myints + sizeof(myints) / sizeof(int) );
    t1    = 10;
    x0    = 0;
    int myints1[] = {x0};
	std::vector<int> x(myints1, myints1 + sizeof(myints1) / sizeof(int) );
    dx0   = 0;
    int myints2[] = {dx0};
	std::vector<int> dx(myints2, myints2 + sizeof(myints2) / sizeof(int) );
    h     = (t1 - t0) / float(N);
    
    for(int i = 1; i<=N; i++) {
        auto diff = RK4(t[i-1],x[i-1],dx[i-1],h);
        x.push_back(  x[i-1] + diff.first  );
        dx.push_back(dx[i-1] + diff.second );
        t.push_back( t[i-1]  + h);
    }
    cout << x[N];
    return 0;
}


which was me trying to replicate the example shown in the vector documentation (http://www.cplusplus.com/reference/vector/vector/vector/). It too failed. I am really trying to fix my own problem, but it seems like I'm missing something as when I asked this at StackOverflow my question was rejected (http://stackoverflow.com/questions/38300582/how-to-add-elements-to-a-vector-in-a-c-loop).
Last edited on
Thomas1965 (4571)
In your first code the error is obvious.
1
2
 
std::vector<int> t;
t[0] = t0;
The vector is empty and you cannot store anything in it. Same with vector x and dx.

Two ways to fix it:
std::vector<int> t(1);
OR
std::vector<int> t;
t.push_back(t0);

fusion809 (3)
Thanks. Although I think there's something a little more subtly wrong with my code. See here is it now:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
 
#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <cmath>
#include <utility>
#include <unistd.h>

using namespace std;

double dx2(int t, int x, int dx)
{
	return (-9.8*cos(x));
}

std::pair<double, double> RK4(float t, float x, float dx, float h)
{
	// k values are diffs in dy/dx.
	// l values are diffs in x.
	double k1, k2, k3, k4, l1, l2, l3, l4, diff1, diff2;
	k1 = h*dx2(t,x,dx);              
	l1 = h*k1;			              
	k2 = h*dx2(t+h/2,x+l1/2,dx+k1/2);
	l2 = h*k2;
	k3 = h*dx2(t+h/2,x+l2/2,dx+k2/2);
	l3 = h*k3;
	k4 = h*dx2(t+h,x+l3,dx+k3);
	l4 = h*k4;
    diff1 = (l1+2*l2+2*l3+l4)/float(6); // diff in x. 
    diff2 = (k1+2*k2+2*k3+k4)/float(6); // diff in y.
    return {diff1, diff2};
}

int main()
{
    double t0, t1, x0, dx0, h, N;
    N     = 100;
    t0    = 0;
    t1    = 10;
    x0    = 0;
    dx0   = 0;
    h     = (t1 - t0) / float(N);
    std::vector<double> t;
    t.push_back(t0);
    std::vector<double> x;
    x.push_back(x0);
    std::vector<double> dx;
    dx.push_back(dx0);
    
    ofstream myfile;
    myfile.open("example.txt");

    for(int i = 1; i<=N; i++) {
        auto diff = RK4(t[i-1],x[i-1],dx[i-1],h);
        t.push_back(  t[i-1] + h);
        x.push_back(  x[i-1] + diff.first  );
        dx.push_back(dx[i-1] + diff.second );
        myfile << diff.first << "\ndiff.first is ";
        myfile << diff.second << "; diff.second is ";
        usleep(1000);
    }
    myfile.close();
    return 0;
}


as you can see I am feeding the x and dx values to my example.txt file which reads:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
 
-0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.098
diff.first is -0.98; diff.second is -0.060458
diff.first is -0.60458; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.0529496
diff.first is -0.529496; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is 0.0251604
diff.first is 0.251604; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is -0.00608362
diff.first is -0.0608362; diff.second is


now what the problem is that as you can see both diff elements are alternating and they should not be. Any ideas why this is?
TheIdeasMan (6817)
Hi,

Looks like a good job for the debugger :+) See where in the code the results differ from what is expected.

Make sure you don't have any divide by zero errors anywhere :+) Avoid using float, prefer double everywhere.

Are your initial values close enough so that they converge? h is 10.0, that might put 2 points on different parts of the curve, giving wild results. On -9.8*cos(x) that would be 10.0 / 2 * pi repetitions of the function.

Just a style issue:

1
2
3
4
5
6
 
double N = 100.0  // put comments here if necessary
    double t0 = 0.0;    // for example valid ranges
    double  t1 = 10.0;
    double x0 = 0.0;
    double dx0 = 0.0;
    double  h = (t1 - t0) / N;  // removed the float cast


This way we guarantee variables are initialised with sensible values.

Hope this helps :+)
fusion809 (3)
Na 10 is the endpoint (t1), 0 is the startpoint, h = (t1-t0)/N; so h = 10/1000 = 0.01. I know that solving this equation, same parameters, with MATLAB goes without a hitch.
TheIdeasMan (6817)
Ok, the debugger then :+)

Naughty Albatross (1283)
It does not look like your program has segfaults in your second solution.

As far as I can see, the same numbers are repeated for a certain number of time. Is it your problem? Is it your goal that you want your values to change constantly during your program run-time?
Topic archived. No new replies allowed.