c - 简谐振子中的改进欧拉法

Question

我已经使用改进的欧拉方法编写了一个 C 代码，以定期确定振荡器的位置、速度和能量。然而，我遇到了一个问题，即振荡器的能量正在减少，尽管没有耗散项。我认为这与我更新位置和速度变量的方式特别相关，并希望就此事获得您的帮助。我的代码如下:

//Compilation and run
//gcc oscillatorimprovedEuler.c -lm -o oscillatorimprovedEuler && ./oscillatorimprovedEuler
#include <stdio.h> 
#include <math.h>  

// The global constans are defined in the following way (having the constant value througout the program
#define m 1.0                 // kg
#define k 1.0                 // kg/sec^2
#define h 0.1                // sec This is the time step
#define N 201 // Number of time steps

int main(void)
{
    // We avoid using arrays this time
    double x = 0, xint = 0;
    double v = 5, vint = 0; // Just like the previous case
    double t = 0;
    double E = (m * v * v + k * x * x) / 2.0; // This is the energy in units of Joules

    FILE *fp = fopen("oscillatorimprovedEuler.dat", "w+");
    int i = 0;
    for(i = 0; i < N ; i++)
    {
        fprintf(fp, "%f \t %f \t %f \t %f \n", x, v, E, t);
        xint = x + (h) * v;
        vint = v - (h) * k * x / m;
        v = v - (h) * ((k * x / m) + (k * xint / m)) / 2.0;
        x = x + (h) * (v + vint) / 2.0;
        E = (m * v * v + k * x * x) / 2.0;
        t += h;
    }
    fclose(fp);
    return 0;
}

可能有一点我遗漏了，如果你能指出来，我将不胜感激。感谢您的帮助。

Answer 1

所以我在 math.stackexchange 的帮助下发现问题与位置和速度的更新时间早于它们应该更新的时间有关，并且需要更多的中间变量。现在的工作代码如下:

//Compilation and run
//gcc oscillatorimprovedEuler.c -lm -o oscillatorimprovedEuler && ./oscillatorimprovedEuler
#include <stdio.h> 
#include <math.h>  

// The global constans are defined in the following way (having the constant value througout the program
#define m 1.0                // kg
#define k 1.0                // kg/sec^2
#define h 0.1                // sec This is the time step
#define N 200                // Number of time steps

int main(void)
{
    // We need to define this many variables to avoid early updating the position and velocity
    double x = 0.0, xpre = 0, xcor = 0;
    double v = 5.0, vpre = 0, vcor = 0; // Just like the previous case
    double t = 0;
    double E = (m * v * v + k * x * x) / 2.0; // This is the energy in units of Joules

    FILE *fp = fopen("oscillatorimprovedEuler.dat", "w+");
    int i = 0;
    for(i = 0; i < N ; i++)
    {
        if (i == 0)
        {
            fprintf(fp, "%f \t %f \t %f \t %f \n", x, v, E, t);
        }
        xpre = x + (h) * v;
        vpre = v - (h) * k * x / m;
        vcor = v - (h) * ((k * x / m) + (k * xpre / m)) / 2.0;
        xcor = x + (h) * (v + vpre) / 2.0;
        E = (m * vcor * vcor + k * xcor * xcor) / 2.0;
        t += h;
        fprintf(fp, "%f \t %f \t %f \t %f \n", xcor, vcor, E, t);
        x = xcor, v = vcor;
    }
    fclose(fp);
    return 0;
}