0

Here is the MATLAB/FreeMat code I got to solve an ODE numerically using the backward Euler method. However, the results are inconsistent with my textbook results, and sometimes even ridiculously inconsistent. What is wrong with the code?

function [x,y] = backEuler(f,xinit,yinit,xfinal,h)

    %f - this is your y prime
    %xinit - initial X
    %yinit - initial Y
    %xfinal - final X
    %h - step size

    n = (xfinal-xinit)/h; %Calculate steps

    %Inititialize arrays...
    %The first elements take xinit and yinit corespondigly, the rest fill with 0s.
    x = [xinit zeros(1,n)];
    y = [yinit zeros(1,n)];

    %Numeric routine
    for i = 1:n
        x(i+1) = x(i)+h;
        ynew = y(i)+h*(f(x(i),y(i)));
        y(i+1) = y(i)+h*f(x(i+1),ynew);
    end
end
4

6 回答 6

5

你的方法是一种新的方法。它既不是后向也不是前向欧拉。:-)

前向欧拉：y1 = y0 + h*f(x0,y0)

后向欧拉solve in y1: y1 - h*f(x1,y1) = y0

你的方法：y1 = y0 +h*f(x0,x0+h*f(x0,y0))

你的方法不是落后的欧拉。

你不解决，你只是用前向欧拉方法y1估计。y1我不想对您的方法进行分析，但我相信即使与前向欧拉相比，它的表现也会很差，因为您f在错误的点评估函数。

这是我能想到的最接近您的方法的方法，也是明确的，应该会产生更好的结果。这是Heun的方法：

y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))

于 2010-05-30T07:29:56.207 回答
2

The only issue I can spot is that the line:

n=(xfinal-xinit)/h

Should be:

n = abs((xfinal-xinit)/h)

To avoid a negative step count. If you are moving in the negative-x direction, make sure to give the function a negative step size.

Your answers probably deviate because of how coarsely you are approximating your answer. To get a semi-accurate result, deltaX has to be very very small and your step size has to be very very very small.

PS. This isn't the "backward Euler method," it is just regular old Euler's method.

If this is homework please tag it so.

于 2010-05-30T01:20:34.343 回答
2

Have a look at numerical recipes, specifically chapter 16, integration of ordinary differential equations. Euler's method is known to have problems:

There are several reasons that Euler’s method is not recommended for practical use, among them, (i) the method is not very accurate when compared to other, fancier, methods run at the equivalent stepsize, and (ii) neither is it very stable

So unless you know your textbook is using Euler's method, I wouldn't expect the results to match. Even if it is, you probably have to use an identical step size to get an identical result.

于 2010-05-30T01:22:36.000 回答
1

除非您真的想通过自己编写的 Euler 方法求解 ODE，否则您应该看看内置的 ODE 求解器。

在旁注中：您不需要x(i)像这样在循环内创建：x(i+1) = x(i)+h;. 相反，您可以简单地编写x = xinit:h:xfinal;. 此外，您可能需要编写n = round(xfinal-xinit)/h);以避免警告。

以下是 MATLAB 实现的求解器。

ode45 基于显式 Runge-Kutta (4,5) 公式，即 Dormand-Prince 对。它是一步求解器——在计算 y(tn) 时，它只需要前一个时间点 y(tn-1) 的解。一般来说，ode45 是大多数问题的第一次尝试应用的最佳函数。

ode23 是 Bogacki 和 Shampine 的显式 Runge-Kutta (2,3) 对的实现。在粗公差和中等刚度的情况下，它可能比 ode45 更有效。与 ode45 一样，ode23 也是一个一步求解器。

ode113 是一个变阶 Adams-Bashforth-Moulton PECE 求解器。在严格的公差和 ODE 文件函数的评估成本特别高的情况下，它可能比 ode45 更有效。ode113 是一个多步求解器——它通常需要几个先前时间点的解来计算当前解。

上述算法旨在解决非刚性系统。如果它们看起来非常慢，请尝试使用下面的刚性求解器之一。

ode15s 是基于数值微分公式 (NDF) 的变阶求解器。可选地，它使用通常效率较低的后向微分公式（BDF，也称为 Gear 方法）。与 ode113 一样，ode15s 是一个多步求解器。当 ode45 失败或效率非常低时尝试 ode15s，并且您怀疑问题是僵硬的，或者在解决微分代数问题时。

ode23s 基于修改后的 2 阶 Rosenbrock 公式。由于它是一步求解器，因此在粗公差下它可能比 ode15s 更有效。它可以解决一些 ode15s 无效的僵硬问题。

ode23t 是使用“自由”插值的梯形规则的实现。如果问题只是中等刚性并且您需要没有数值阻尼的解决方案，请使用此求解器。ode23t 可以解决 DAE。

ode23tb 是 TR-BDF2 的一个实现，它是一个隐式 Runge-Kutta 公式，第一阶段是梯形规则步骤，第二阶段是二阶反向微分公式。通过构造，相同的迭代矩阵用于评估两个阶段。与 ode23s 一样，此求解器在粗公差下可能比 ode15s 更有效。

于 2010-05-30T01:40:45.417 回答
0

代码很好。只是你必须在 for 循环中添加另一个循环。检查一致性水平。

if abs((y(i+1) - ynew)/ynew) > 0.0000000001 
    ynew = y(i+1);
    y(i+1) = y(i)+h*f(x(i+1),ynew);
end

我检查了一个虚拟函数，结果很有希望。

于 2016-05-17T15:38:50.670 回答
0

我认为这段代码可以工作。尝试这个。

for i =1:n
    t(i +1)=t(i )+dt;
    y(i+1)=solve('y(i+1)=y(i)+dt*f(t(i+1),y(i+1)');
    end
于 2016-03-15T04:37:30.827 回答