#include #include #include #include using namespace std; typedef std::vector vec; typedef std::vector Matrix; Matrix Hein(double a, double h, double ya, double dya, double p, double q, std::function fx) { Matrix Y(4); std::function 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 fx) { double h = (b - a) / (n - 1); Matrix Y(n); vec F(n); std::function 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 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; } }