MATLAB代码帮助。向后的Euler方法

StackOverflow https://stackoverflow.com/questions/2937183

题

这里是 MATLAB/Freemat 我要解决的代码 颂 数值使用 向后的欧拉 方法。但是，结果与我的教科书结果不一致，有时甚至不一致。代码怎么了？

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
有帮助吗？

解决方案

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

前锋欧拉： 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 使用正向Euler方法。我不想继续对您的方法进行分析，但我认为即使与前锋欧拉相比，它的确会表现不佳，因为您评估了功能 f 在错误的点。

这是我能想到的最接近您的方法的方法，它也应该给您带来更好的结果。它是 Heun的方法:

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

其他提示

我唯一可以发现的问题是这条线：

n=(xfinal-xinit)/h

应该：

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

为了避免否定的步骤计数。如果您沿负X方向移动，请确保使函数的步长为负。

您的答案可能会偏离，因为您对答案的近似程度很高。要获得半准确的结果，Deltax必须非常小，您的步骤尺寸必须非常非常小。

PS。这不是“向后的欧拉方法”，它只是常规的旧欧拉方法。

如果这是作业，请标记。

看一下 数值食谱, ，特别是第16章，是普通微分方程的集成。已知Euler的方法有问题：

不建议将Euler的方法用于实际使用的原因有几个，（i）与其他相比，该方法在等效步骤中运行的方法不是很准确，并且（ii）也不是很稳定的。

因此，除非您知道您的教科书正在使用Euler的方法，否则我不会期望结果匹配。即使是，您也可能必须使用相同的步长才能获得相同的结果。

除非你 真的 想要通过您自己写的Euler的方法来解决颂歌，您应该看看 内置的ODE求解器.

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

这是MATLAB实施的求解器。

ODE45基于Dormand-prince对的显式runge-kutta（4,5）公式。它是一个单步求解器 - 在计算Y（TN）中，仅需要在即将到达的时间点y（TN-1）的解决方案。通常，ODE45是针对大多数问题的首次尝试应用的最佳功能。

ODE23是一对bogacki和洗发剂的显式runge-kutta（2,3）的实现。在粗耐受性和中等刚度的情况下，它可能比Ode45更有效。像ode45一样，ode23是一个单步求解器。

ODE113是可变订单Adams-Bashforth-Moulton Pece求解器。在严格的公差时，它可能比Ode45更有效，并且何时评估ODE文件功能特别昂贵。 ODE113是一个多步求解器 - 通常需要在上一个时间点的几个时间点来计算当前解决方案。

上述算法旨在求解非固定系统。如果它们似乎过慢，请尝试使用下面的刚性求解器之一。

ODE15S是基于数值分化公式（NDFS）的可变订单求解器。可选地，它使用效率较低的后退分化公式（BDF，也称为Gear's方法）。像ODE113一样，ODE15S是多步求解器。当ODE45失败或非常低效时，请尝试使用ODE15S，并且您怀疑问题很严重，或者在解决差异代数问题时。

ODE23S基于订单2的修改后的Rosenbrock公式。由于它是一个单步求解器，因此在粗公也可以比Ode15更有效。它可以解决ode15s无效的某些严格问题。

ODE23T是使用“自由”插值的梯形规则的实现。如果问题仅适度僵硬，并且您需要一个无数值阻尼的解决方案，请使用该求解器。 ODE23T可以解决DAE。

ODE23TB是TR-BDF2的实现，TR-BDF2是一个隐含的runge-kutta公式，其第一阶段是梯形规则步骤和第二阶段，第二阶段是第二阶的向后分化公式。通过构造，使用相同的迭代矩阵来评估两个阶段。像ODE23S一样，该求解器可能比在粗耐受性下的ODE15更有效。

我认为此代码可以起作用。尝试这个。

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

代码很好。只要您必须在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

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

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top