C++

Forum

 
 
runge-kutta with 2d kinematics

crimsonzero2 (213)
I got to make a runge-kutta program that can solve where an object is at an x position and a y position. Where dx/dt = vx, dy/dt = vy, dvx/dt = 0, and dvy/dt = -g. The initial launch is at 20 m/s with a 30 degrees angle. X and y are suppose to start a 0. The equation for x is x = vx0 * t and y is y = vy0 * t - (1/2) * g * t^2. I am having a problem with how to use the initial condition of vx and vy with the rate of change it has. Here is the code I got so far:

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
102
103
 
#include <iostream>
#include <cmath>

double ddt[4];

void f(double t, double v[4])
{
	ddt[0] = v[2] * t;
	ddt[1] = v[3] * t - (0.5) * 9.81 * pow(t,2);
	ddt[2] = 0 * t;
	ddt[3] = -9.81 * t;
}

int main()
{
	using namespace std;
	void f(double t, double v[4]);
	
	double h, t0, v1_0, v2_0, v3_0, v4_0;
	int i, j;
	const int nos = 10;
	
	h = 0.1;
	t0 = 0.0;
	v1_0 = 0.0;
	v2_0 = 0.0;
	v3_0 = 20 * cos(30);
	v4_0 = 20 * sin(30);
	
	double t[nos + 1];
	double v[4];
	
	t[0] = t0;
	v[0] = v1_0;
	v[1] = v2_0;
	v[2] = v3_0;
	v[3] = v4_0;
	
	double v_n[4];
	
	cout << "x is: " << v[0] << endl;
	cout << "y is: " << v[1] << endl << endl;
	
	for(i = 0; i < nos; i++)
	{
		t[i + 1] = t[i] + h;
		v_n[0] = v[0];
		v_n[1] = v[1];
		v_n[2] = v[2];
		v_n[3] = v[3];
		
		double added[4];
		
		double k1[4], k2[4], k3[4], k4[4];
		double v_npl[4];
		
		f(t[i], v_n);
		for(j = 0; j < 4; j++)
		{
			k1[j] = h * ddt[j];
			added[j] = v_n[j] + (0.5 * k1[j]);
		}
		
		f(t[i] = 0.5 * h, added);
		for(j = 0; j < 4; j++)
		{
			k2[j] = h * ddt[j];
			added[j] = v_n[j] + (0.5 * k2[j]);
		}
		
		f(t[i] + 0.5 * h, added);
		for(j = 0; j < 4; j++)
		{
			k3[j] = h * ddt[j];
			added[j] = v_n[j] + k3[j];
		}
		
		f(t[i] + h, added);
		for(j = 0; j < 4; j++)
		{
			k4[j] = h * ddt[j];
		}
		
		for(j = 0; j < 4; j++)
		{
			v_npl[j] = v_n[j] + ((k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0);
		}
		
		v[0] = v_npl[0];
		v[1] = v_npl[1];
		v[2] = v_npl[2];
		v[3] = v_npl[3];
		
		cout << "x is: " << v[0] << endl;
		cout << "y is: " << v[1] << endl << endl;
	}
	
	return 0;
}

/*
doesn't seem to give "+" to y and y never returns to 0
*/


And here is the output:


x is: 0
y is: 0

x is: 0.0231377
y is: -0.152415

x is: 0.0514171
y is: -0.339338

x is: 0.0848383
y is: -0.562826

x is: 0.123401
y is: -0.825019

x is: 0.167106
y is: -1.12814

x is: 0.215952
y is: -1.47448

x is: 0.26994
y is: -1.86644

x is: 0.32907
y is: -2.30648

x is: 0.393341
y is: -2.79714

x is: 0.462754
y is: -3.34106


The problem I'm having is that y is negative and it won't return to 0. Anyone know what I am doing wrong? Thank you helping me if you do.
crimsonzero2 (213)
It was a simple radian-degree error. Here is the 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
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
102
 
#include <iostream>
#include <cmath>

#define PI 3.1415926

double ddt[4];

void f(double t, double v[4])
{
	ddt[0] = v[2] * t;
	ddt[1] = v[3] * t - (0.5) * 9.81 * pow(t,2);
	ddt[2] = 0;
	ddt[3] = -9.81;
}

int main()
{
	using namespace std;
	void f(double t, double v[4]);
	
	double h, t0, v1_0, v2_0, v3_0, v4_0;
	int i, j;
	const int nos = 10;
	
	h = 0.1;
	t0 = 0.0;
	v1_0 = 0.0;
	v2_0 = 0.0;
	v3_0 = 20 * cos(30 * PI/180);
	v4_0 = 20 * sin(30 * PI/180);
	
	double t[nos + 1];
	double v[4];
	
	t[0] = t0;
	v[0] = v1_0;
	v[1] = v2_0;
	v[2] = v3_0;
	v[3] = v4_0;
	
	double v_n[4];
	
	cout << "x is: " << v[0] << endl;
	cout << "y is: " << v[1] << endl << endl;
	
	for(i = 0; i < nos; i++)
	{
		t[i + 1] = t[i] + h;
		v_n[0] = v[0];
		v_n[1] = v[1];
		v_n[2] = v[2];
		v_n[3] = v[3];
		
		double added[4];
		
		double k1[4], k2[4], k3[4], k4[4];
		double v_npl[4];
		
		f(t[i], v_n);
		for(j = 0; j < 4; j++)
		{
			k1[j] = h * ddt[j];
			added[j] = v_n[j] + (0.5 * k1[j]);
		}
		
		f(t[i] = 0.5 * h, added);
		for(j = 0; j < 4; j++)
		{
			k2[j] = h * ddt[j];
			added[j] = v_n[j] + (0.5 * k2[j]);
		}
		
		f(t[i] + 0.5 * h, added);
		for(j = 0; j < 4; j++)
		{
			k3[j] = h * ddt[j];
			added[j] = v_n[j] + k3[j];
		}
		
		f(t[i] + h, added);
		for(j = 0; j < 4; j++)
		{
			k4[j] = h * ddt[j];
		}
		
		for(j = 0; j < 4; j++)
		{
			v_npl[j] = v_n[j] + ((k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0);
		}
		
		v[0] = v_npl[0];
		v[1] = v_npl[1];
		v[2] = v_npl[2];
		v[3] = v_npl[3];
		
		cout << "x is: " << v[0] << endl;
		cout << "y is: " << v[1] << endl << endl;
	}
	
	return 0;
}
Topic archived. No new replies allowed.