Решить уравнение при помощи методов Эйлера и Рунге Кутта

[Gladiator1995]

Здравсвуйте, дамы и господа! В институте дали такое задание:
Вот условие задачи
Решить методом Эйлера и Рунге-Кутта уравнение, описывающее смещение x электрона в атоме под действием электрического поля монохроматической световой волны

$$\ddot{x} + \frac{1}{\tau}\dot{x} + {\omega }^{2}x + \gamma {x}^{2} = eE_{0}\cos (pt)$$

\omega - Круговая частота собственных колебаний электрона в атоме
\tau - Их время затухания
\gamma - Коэффициент нелинейности возвращающей силы, действующей на электрон со стороны ядра
e, m - заряд и масса электрона
E_{0} - амплитудное значение напряженности электрического поля световой волны
p - круговая частота волны

Написал программу, касающуюся метода Эйлера, но ее фазовый портрет отвратителен! Помогите найти ошибку

Используется синтаксис C

<img src='https://dxdy-04.korotkov.co.uk/f/7/7/a/77a90a0f3ce59a721411fa13cce0768782.png'  class="latex" alt="#ifndef FUNCTION
#define FUNCTION

double F1(double E0, double p, double t, double tau, double omega, double gamma, double x, double y);

double F2(double E0, double p, double t, double tau, double omega, double gamma, double x, double y);

#endif" title="#ifndef FUNCTION
#define FUNCTION

double F1(double E0, double p, double t, double tau, double omega, double gamma, double x, double y);

double F2(double E0, double p, double t, double tau, double omega, double gamma, double x, double y);

#endif" />

<img src='https://dxdy-04.korotkov.co.uk/f/7/b/0/7b0132310eef7b9bb7d934e13ec3269f82.png'  class="latex" alt="#include "Function.h"
#include <stdio.h>
#include <math.h>

#define q 1,6E-19
#define m 9,1E-27

double F1(double E0, double p, double t, double tau, double omega, double gamma, double x, double y)
{
        // t - время
        // E0 - амплитудное значение напряженности электрического поля световой волны
        // р - круговая частота волны
        // tau - время затухания свободных колебаний
        // omega - круговая частота свободных колебаний электрона в атоме
        // gamma - коэффициент нелинейности возвращающей силы
        return (q * E0 * cos(p * t)/m - y/(tau) - (omega * omega * x) - gamma * x * x); // F1(x, y, t)
}

double F2(double E0, double p, double t, double tau, double omega, double gamma, double x, double y)
{
        // t - время
        // E0 - амплитудное значение напряженности электрического поля световой волны
        // р - круговая частота волны
        // tau - время затухания свободных колебаний
        // omega - круговая частота свободных колебаний электрона в атоме
        // gamma - коэффициент нелинейности возвращающей силы
        return y; // F2(x, y, t)
}" title="#include "Function.h"
#include <stdio.h>
#include <math.h>

#define q 1,6E-19
#define m 9,1E-27

double F1(double E0, double p, double t, double tau, double omega, double gamma, double x, double y)
{
        // t - время
        // E0 - амплитудное значение напряженности электрического поля световой волны
        // р - круговая частота волны
        // tau - время затухания свободных колебаний
        // omega - круговая частота свободных колебаний электрона в атоме
        // gamma - коэффициент нелинейности возвращающей силы
        return (q * E0 * cos(p * t)/m - y/(tau) - (omega * omega * x) - gamma * x * x); // F1(x, y, t)
}

double F2(double E0, double p, double t, double tau, double omega, double gamma, double x, double y)
{
        // t - время
        // E0 - амплитудное значение напряженности электрического поля световой волны
        // р - круговая частота волны
        // tau - время затухания свободных колебаний
        // omega - круговая частота свободных колебаний электрона в атоме
        // gamma - коэффициент нелинейности возвращающей силы
        return y; // F2(x, y, t)
}" />

<img src='https://dxdy-01.korotkov.co.uk/f/0/8/2/082bda10334f645d7c1390e23065a2af82.png'  class="latex" alt="#ifndef EULER
#define EULER
#include <stdio.h>

void Euler(FILE *ef, double E0, double p, double t, double tau, double omega, double gamma, const double a, const double b, const double eps);

#endif" title="#ifndef EULER
#define EULER
#include <stdio.h>

void Euler(FILE *ef, double E0, double p, double t, double tau, double omega, double gamma, const double a, const double b, const double eps);

#endif" />

<img src='https://dxdy-03.korotkov.co.uk/f/e/7/2/e72121ebdb2676ca5b8fa37ab1fcff8782.png'  class="latex" alt="#include "Function.h"
#include "Euler.h"
#include <math.h>

void Euler(FILE *ef, double E0, double p, double t, double tau, double omega, double gamma, const double a, const double b, const double eps)
{
        double x, x2, y, y2, dt;
        x = a;
        y = b;
        dt = 0.01;
        fprintf(ef, "Initial data x(0) = a = %f; y(0) = b = %f; Epsilon = eps = %f;\n", x, y, eps);
        fprintf(ef, "E0 = %f;", E0);
        fprintf(ef, "\n");
        fprintf(ef, "p = %f;", p);
        fprintf(ef, "\n");
        fprintf(ef, "tau = %f;", tau);
        fprintf(ef, "\n");
        fprintf(ef, "omega = %f;", omega);
        fprintf(ef, "\n");
        fprintf(ef, "gamma = %f;", gamma);
        for(int i = 0;i < 10000; ++i)
                {       
                        while (fabs((F1(E0, p, t, tau, omega, gamma, x, y) * dt / 2)) > eps || fabs((F2(E0, p, t, tau, omega, gamma, x, y) * dt / 2)) > eps)    dt/=2.0;
                        x2 = x + F1(E0, p, t, tau, omega, gamma, x, y) * dt; // ???????? x ? y ? ???? ?????? ?????? ???? ?? ?? ...
                        y2 = y + F2(E0, p, t, tau, omega, gamma, x, y) * dt; // ... ??? ? ? ???? ??????
                        x = x2;
                        y = y2;
                        t += dt;
                        fprintf(ef, "%f   %f   %f\n", t , x2, y2);
                        dt = 0.01;
                }
}" title="#include "Function.h"
#include "Euler.h"
#include <math.h>

void Euler(FILE *ef, double E0, double p, double t, double tau, double omega, double gamma, const double a, const double b, const double eps)
{
        double x, x2, y, y2, dt;
        x = a;
        y = b;
        dt = 0.01;
        fprintf(ef, "Initial data x(0) = a = %f; y(0) = b = %f; Epsilon = eps = %f;\n", x, y, eps);
        fprintf(ef, "E0 = %f;", E0);
        fprintf(ef, "\n");
        fprintf(ef, "p = %f;", p);
        fprintf(ef, "\n");
        fprintf(ef, "tau = %f;", tau);
        fprintf(ef, "\n");
        fprintf(ef, "omega = %f;", omega);
        fprintf(ef, "\n");
        fprintf(ef, "gamma = %f;", gamma);
        for(int i = 0;i < 10000; ++i)
                {       
                        while (fabs((F1(E0, p, t, tau, omega, gamma, x, y) * dt / 2)) > eps || fabs((F2(E0, p, t, tau, omega, gamma, x, y) * dt / 2)) > eps)    dt/=2.0;
                        x2 = x + F1(E0, p, t, tau, omega, gamma, x, y) * dt; // ???????? x ? y ? ???? ?????? ?????? ???? ?? ?? ...
                        y2 = y + F2(E0, p, t, tau, omega, gamma, x, y) * dt; // ... ??? ? ? ???? ??????
                        x = x2;
                        y = y2;
                        t += dt;
                        fprintf(ef, "%f   %f   %f\n", t , x2, y2);
                        dt = 0.01;
                }
}" />

<img src='https://dxdy-02.korotkov.co.uk/f/d/4/1/d41d8cd98f00b204e9800998ecf8427e82.png'  class="latex" alt="#include "Euler.h"
#include <stdio.h>

#define q 1,6E-19
#define m 9,1E -27

int main()
{
        double x, y, E0, p, t, tau, omega, gamma, eps;
        // t - âðåìÿ
        // E0 - àìïëèòóäíîå çíà÷åíèå íàïðÿæåííîñòè ýëåêòðè÷åñêîãî ïîëÿ ñâåòîâîé âîëíû
        // ð - êðóãîâàÿ ÷àñòîòà âîëíû
        // tau - âðåìÿ çàòóõàíèÿ ñâîáîäíûõ êîëåáàíèé
        // omega - êðóãîâàÿ ÷àñòîòà ñâîáîäíûõ êîëåáàíèé ýëåêòðîíà â àòîìå
        // gamma - êîýôôèöèåíò íåëèíåéíîñòè âîçâðàùàþùåé ñèëû
        // epsilon = eps - 
        printf("Enter epsilon\n");
        scanf("%lf", &eps);
        printf("Enter x, y\n");
        scanf("%lf %lf", &x, &y);
        printf("Enter E0:\n");
        scanf("%lf", &E0);
        printf("Enter p:\n");
        scanf("%lf", &p);
        printf("Enter t:\n");
        scanf("%lf", &t);
        printf("Enter tau:\n");
        scanf("%lf", &tau);
        printf("Enter omega:\n");
        scanf("%lf", &omega);
        printf("Enter gamma:\n");
        scanf("%lf", &gamma);
        printf("\n");
        FILE *ef;
        ef = fopen("Euler.txt", "w");
        Euler(ef, E0, p, t, tau, omega, gamma, x, y, eps);
        fclose(ef);
        return 0;
}" title="#include "Euler.h"
#include <stdio.h>

#define q 1,6E-19
#define m 9,1E -27

int main()
{
        double x, y, E0, p, t, tau, omega, gamma, eps;
        // t - âðåìÿ
        // E0 - àìïëèòóäíîå çíà÷åíèå íàïðÿæåííîñòè ýëåêòðè÷åñêîãî ïîëÿ ñâåòîâîé âîëíû
        // ð - êðóãîâàÿ ÷àñòîòà âîëíû
        // tau - âðåìÿ çàòóõàíèÿ ñâîáîäíûõ êîëåáàíèé
        // omega - êðóãîâàÿ ÷àñòîòà ñâîáîäíûõ êîëåáàíèé ýëåêòðîíà â àòîìå
        // gamma - êîýôôèöèåíò íåëèíåéíîñòè âîçâðàùàþùåé ñèëû
        // epsilon = eps - 
        printf("Enter epsilon\n");
        scanf("%lf", &eps);
        printf("Enter x, y\n");
        scanf("%lf %lf", &x, &y);
        printf("Enter E0:\n");
        scanf("%lf", &E0);
        printf("Enter p:\n");
        scanf("%lf", &p);
        printf("Enter t:\n");
        scanf("%lf", &t);
        printf("Enter tau:\n");
        scanf("%lf", &tau);
        printf("Enter omega:\n");
        scanf("%lf", &omega);
        printf("Enter gamma:\n");
        scanf("%lf", &gamma);
        printf("\n");
        FILE *ef;
        ef = fopen("Euler.txt", "w");
        Euler(ef, E0, p, t, tau, omega, gamma, x, y, eps);
        fclose(ef);
        return 0;
}" />

[Pphantom]

Тема перемещена из форума «Программирование» в форум «Карантин»
по следующим причинам:

- неправильно набраны формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы);
- неправильно оформлен код.

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.