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;
}