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>
//function 1 (u1_0 uses this function as u1_1 uses the other function and t1 is time)
double f1(double t1, double u1_0, double u1_1)
{
double dudt1 = (3 * u1_0) + (2 * u1_1) - (2 * (t1 * t1) + 1) * exp(2 * t1);
return dudt1;
}
//function 2 (u2_0 uses the other function as u2_1 uses this function and t2 is time)
double f2(double t2, double u2_0, double u2_1)
{
double dudt2 = (4 * u2_0) + u2_1 - ((t2 * t2) + (2 * t2) - 4) * exp(2 * t2);
return dudt2;
}
int main ()
{
using namespace std;
//variables
double h, t0, u1_0, u2_0;
//loop variable and number of steps
int i, nos;
nos = 10;
h = 0.1;
t0 = 0.0;
u1_0 = 1.0;
u2_0 = 1.0;
//array variables
double t[nos + 1];
double u[2][nos + 1];
double u_n[2][1];
t[0] = t0;
u[0][0] = u1_0; //u[0][0] = 1.0
u[1][0] = u2_0; //u[1][0] = 1.0
double u_npl[2][1];
//Runge-Kutta variables for both functions above
double k1[2], k2[2], k3[2], k4[2];
cout << "This is the RK4 Method: " << endl << endl;
for(i = 0; i <= nos; i++)
{
t[i + 1] = t[i] + h;
u_n[0][0] = u[0][i];
u_n[1][0] = u[1][i];
k1[0] = h * f1(t[i], u_n[0][0], u_n[1][0]);
k1[1] = h * f2(t[i], u_n[0][0], u_n[1][0]);
//not sure if this is how you would get k2 with three variables
k2[0] = h * f1(t[i] + 0.5 * h, u_n[0][0] + 0.5 * k1[0], u_n[1][0] + 0.5 * k1[0]);
k2[1] = h * f2(t[i] + 0.5 * h, u_n[0][0] + 0.5 * k1[1], u_n[1][0] + 0.5 * k1[1]);
//not sure if this is how you would get k3 with three variables
k3[0] = h * f1(t[i] + 0.5 * h, u_n[0][0] + 0.5 * k2[0], u_n[1][0] + 0.5 * k2[0]);
k3[1] = h * f2(t[i] + 0.5 * h, u_n[0][0] + 0.5 * k2[1], u_n[1][0] + 0.5 * k2[1]);
k4[0] = h * f1(t[i] + h, u_n[0][0] + k3[0], u_n[1][0] + k3[0]);
k4[1] = h * f2(t[i] + h, u_n[0][0] + k3[1], u_n[1][0] + k3[1]);
u_npl[0][0] = u_n[0][0] + ((k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6);
u_npl[1][0] = u_n[1][0] + ((k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6);
u[0][i + 1] = u_npl[0][0];
u[1][i + 1] = u_npl[1][0];
cout << "value of u[0][" << i << "] is: " << u[0][i] << endl;
cout << "value of u[1][" << i << "] is: " << u[1][i] << endl << endl;
}
cout << "This is the actual result: " << endl << endl;
u[0][0] = u1_0;
u[1][0] = u2_0;
t[0] = t0;
for (i = 0; i <= nos; i++)
{
t[i + 1] = t[i] + h;
u_npl[0][0] = (0.333333 * exp( 5 * t[i])) - (0.333333 * exp(-t[i])) + exp(2 * t[i]);
u_npl[1][0] = (0.333333 * exp(5 * t[i])) + (0.666667 * exp(-t[i])) + ((t[i] * t[i]) * exp( 2 * t[i]));
u[0][i + 1] = u_npl[0][0];
u[1][i + 1] = u_npl[1][0];
cout << "value of u[0][" << i << "] is: " << u[0][i] << endl;
cout << "value of u[1][" << i << "] is: " << u[1][i] << endl << endl;
}
return 0;
}
|
