我的欧拉方法实现是否正确?

Is my euler's method implementation correct?

我正在编写代码来求解已经使用欧拉方法给出的微分方程。得到结果后,我将它们绘制在图表上。但是我得到了不正确的值。

这是为了得到t运行sfer函数的反馈响应。我已经尝试以其他方式编写 euler 的方法(结果相同)。如果我的微分方程写得正确,我检查了几次。我想知道问题出在这里还是其他地方。我 运行 不知道这段代码有什么问题。

public static double inputSignal(double t, double A, double f, int typeOfSignal) {

        switch (typeOfSignal) {
            case 1:
                return A*cos(2*Math.PI*f*t);
            case 2:
                return A*(2/Math.PI)*asin(sin(2*Math.PI*f*t));
            case 3:
                return A*signum(sin(2*Math.PI*f*t));
            default:
                break;
        }

        return 0;
    }

public static double equation(double x, double xDerivative, double u, double uDerivative, double a0, double a1, double b0, double b1, double b2) {

        return (a1*uDerivative + a0*u - (a0+b0)*x - (a1+b1)*xDerivative)/b2; 
    }

 public static void eulerMethod(double h, double A, double f, double a0, double a1, double b0, double b1, double b2, int typeOfSignal) {

        double xn = 0, yn = 0, tn = 0, tn1; //x' = y; y' = f(x, y, u)
        double uDerivative;

        for (int i = 0; i < 10000; i++) {
            tn1 = tn + h;
            uDerivative= (inputSignal(tn1+h, A, f, typeOfSignal)-inputSignal(tn1, A, f, typeOfSignal))/h;
            xn = xn + h*equation(xn, yn, pobudzenie(tn1, A, f, typeOfSignal), uDerivative, a0, a1, b0, b1, b2);
            yn = yn + h*xn;
            System.out.println("Step" + i);
            System.out.println("Y(" + tn1 + ") " + " = " + yn);
            tn = tn1;
            results[i] = yn;
            inputGraph[i] = inputSignal(tn, A, f, typeOfSignal);
        }

    }

微分方程为:b2*x'' = a1*u' + a0*u - (a0+b0)*x - (a1+b1)*x'

我想 euler2(0.0025, 1, 0.01, 1, 5, 1, 2, 1, 1) 得到的是 this, but what I get from my code is this.

如您所见,我的函数有不想要的波动,它的振幅低于预期(它是 0.2,而它应该是 0.7)。

最有可能的错误

你在你的一阶系统中混淆了一些东西。你有

x'' = f(t,x,x',u,u')

和设置 y=x' 这应该导致第一顺序系统

x' = y
y' = f(t,x,y,u,u')

但是,您在计算 xy 的更新时交换了右侧,有效地将等式更改为 x''=f(t,x',x,u,u') 并显示 x'

在您的 semi-explicit Euler 实现中,正确的方法实现是

yn = yn + h*equation(xn, yn, inputSignal(tn1, A, f, typeOfSignal), uDerivative, a0, a1, b0, b1, b2);
xn = xn + h*yn;

在严格的显式 Euler 实现中,您不会在 xn 更新中使用新计算的 yn,而是使用旧值。需要涉及一些临时缓存变量。


改进的一阶系统

您可以通过设置

获得更直接的 first-order 系统
y=b2*x'+(a1+b1)*x-a1*u

这样

x' = ( y - (a1+b1)*x+a1*u )/b2;
y' = a0*u-(a0+b0)*x

从 first-order 系统中删除 u 的所有衍生物。

Visual ODE 求解器测试

如果您为任一系统生成具有不同步长的图,则视觉收敛表明实现正确。那么控制图的任何差异都意味着传递给此处的求解器和传递给控制求解器的 ODE 可能存在差异。