include iostream include cmath include functional include vector using

 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
#include <iostream>
#include <cmath>
#include <functional>
#include <vector>
using namespace std;
typedef std::vector<double> vec;
typedef std::vector<vec> Matrix;
Matrix Hein(double a, double h, double ya, double dya, double p, double q, std::function<double(double)> fx)
{
Matrix Y(4);
std::function<vec(double, vec)> f = [&](double xx, vec yy) -> vec{
vec temp ={ yy[1], fx(xx) - p * yy[1] - q * yy[0]};
return temp;
};
const int m = 10;
vec curPt = { ya, dya };
h /= m;
for (int i = 0; i <= 3 * m; ++i)
{
double x = a + i * h;
if (i % m == 0) Y[i / m] = curPt;
vec fi = f(x, curPt);
vec newPt = { curPt[0] + h * fi[0], curPt[1] + h * fi[1] };
vec fii = f(x + h, newPt);
newPt = {
curPt[0] + h / 2 * (fi[0] + fii[0]),
curPt[1] + h / 2 * (fi[1] + fii[1]),
};
curPt = newPt;
}
return Y;
}
Matrix Solve(double a, double b, int n, double ya, double dya, double p, double q, std::function<double(double)> fx)
{
double h = (b - a) / (n - 1);
Matrix Y(n);
vec F(n);
std::function<double(double, double, double)> f = [&](double x, double y, double dy) -> double{
return fx(x) - p * dy - q * y;
};
Matrix StartPts = Hein(a, h, ya, dya, p, q, fx);
for (int i = 0; i < 4; ++i)
{
double x = a + i * h;
Y[i] = StartPts[i];
F[i] = f(x, Y[i][0], Y[i][1]);
}
for (int i = 4; i < n; ++i)
{
double x = a + i * h;
Y[i] =
{
Y[i - 1][0] + h * Y[i - 1][1] + h * h / 360 * (323 * F[i - 1] - 264 * F[i - 2] + 159 * F[i - 3] - 38 * F[i - 4]),
Y[i - 1][1] + h / 24 * (55 * F[i - 1] - 59 * F[i - 2] + 37 * F[i - 3] - 9 * F[i - 4])
};
F[i] = f(x, Y[i][0], Y[i][1]);
}
return Y;
}
int main()
{
double p = -3, q = 2, A = 4, B = -1, a = 0.15, b = 0.85, h = 0.05;
int n = std::roundf((b - a) / h ) + 1;
std::function<double(double)> fx = [&](double x) -> double
{
return x*x - 7*x;
};
Matrix temp = Solve(a, b, n, A, 0.13585390768417982, p, q, fx);
for (int i = 0; i < n; ++i)
{
cout << "y(" << a+i*h << ")" << " = " << temp[i][0] << endl;
}
}