MATLAB 代码帮助。后向欧拉法

发布于 2024-09-03 15:11:22 字数 1005 浏览 9 评论 0原文

这是 MATLAB/FreeMat 代码我必须解决 ODE< /a> 使用向后欧拉方法进行数值计算。然而,结果与我的课本结果不一致,有时甚至不一致得可笑。代码有什么问题吗?

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

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

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(6

向日葵 2024-09-10 15:11:22

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

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

后向欧拉求解 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)))

Your method is a method of a new kind. It is neither backward nor forward Euler. :-)

Forward Euler: y1 = y0 + h*f(x0,y0)

Backward Euler solve in y1: y1 - h*f(x1,y1) = y0

Your method: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Your method is not backward Euler.

You don't solve in y1, you just estimate y1 with the forward Euler method. I don't want to pursue the analysis of your method, but I believe it will behave poorly indeed, even compared with forward Euler, since you evaluate the function f at the wrong point.

Here is the closest method to your method that I can think of, explicit as well, which should give much better results. It's Heun's Method:

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

安稳善良 2024-09-10 15:11:22

我能发现的唯一问题是这一行:

n=(xfinal-xinit)/h

应该是:

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

避免负步数。如果您向负 x 方向移动,请确保为函数指定负步长。

您的答案可能会有所偏差,因为您对答案的近似程度有多粗略。为了获得半准确的结果,deltaX 必须非常非常小,并且步长必须非常非常小。

附言。这不是“向后欧拉方法”,它只是常规的旧欧拉方法。

如果这是家庭作业,请标记它。

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.

万劫不复 2024-09-10 15:11:22

看看数值食谱,特别是第16章,常微分方程的积分。众所周知,欧拉方法存在问题:

不推荐欧拉方法实际使用有几个原因,其中,(i) 与以等效步长运行的其他更高级的方法相比,该方法不是很准确,(ii) 也不是很准确非常稳定

稳定除非您知道您的教科书使用的是欧拉方法,否则我不希望结果匹配。即使是这样,您可能也必须使用相同的步长才能获得相同的结果。

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.

音盲 2024-09-10 15:11:22

除非您确实想通过自己编写的欧拉方法求解 ODE,否则您应该看看 内置 ODE 求解器

旁注:您不需要在循环内创建 x(i),如下所示:x(i+1) = x(i)+h;。相反,您可以简单地编写 x = xinit:h:xfinal;。另外,您可能需要编写 n = round(xfinal-xinit)/h); 以避免警告。


以下是 MATLAB 实现的求解器。

ode45 基于显式
龙格-库塔(4,5)公式,
多蒙德王子对。这是一步到位
求解器 – 在计算 y(tn) 时,需要
只有立即解决
之前的时间点 y(tn-1)。在
一般来说,ode45 是最好的函数
对于大多数人来说作为第一次尝试申请
问题。

ode23 是一个实现
显式龙格-库塔 (2,3) 对
博加基和香平。可能会更多
粗加工时比 ode45 高效
公差和存在
刚度适中。像 ode45、ode23
是一个一步求解器。

ode113 是一个可变阶数
Adams-Bashforth-Moulton PECE 求解器。
它可能比 ode45 更高效
严格的公差和当 ODE
文件功能特别
评估昂贵。 ode113 是一个
多步求解器——通常需要
前面几处的解决方案
计算电流的时间点
解决方案。

上述算法的目的是
解决非刚性系统。如果他们出现
如果速度太慢,请尝试使用其中之一
下面的刚性求解器。

ode15s 是一个变阶求解器
基于数值微分
公式(NDF)。可选地,它使用
后向微分公式
(BDF,也称为 Gear 方法)
通常效率较低。喜欢
ode113、ode15s 是多步求解器。
当 ode45 失败时尝试 ode15s,或者
效率非常低,你怀疑
问题很僵化,或者在解决时
微分代数问题。

ode23s是基于修改的
Rosenbrock 2 阶公式。因为
它是一个一步求解器,它可能是
在原油条件下比 ode15 更高效
公差。它可以解决某些类型的
ode15s 无法解决的刚性问题
有效的。

ode23t 是一个实现
使用“自由”的梯形规则
插值。如果满足以下条件,请使用此求解器
问题只是中等硬度并且
你需要一个没有数值的解决方案
减震。 ode23t 可以解决 DAE。

ode23tb 是一个实现
TR-BDF2,隐式龙格库塔
公式的第一阶段是
梯形法则步骤和第二步
落后的阶段
二阶微分公式。
通过构造,同样的迭代
矩阵用于评估两者
阶段。与 ode23s 一样,该求解器可以
在原油条件下比 ode15 更高效
公差。

Unless you really want to solve an ODE via Euler's method that you've written by yourself you should have a look at built-in ODE solvers.

On a sidenote: you don't need to create x(i) inside the loop like this: x(i+1) = x(i)+h;. Instead, you can simply write x = xinit:h:xfinal;. Also, you may want to write n = round(xfinal-xinit)/h); to avoid warnings.


Here are the solvers implemented by MATLAB.

ode45 is based on an explicit
Runge-Kutta (4,5) formula, the
Dormand-Prince pair. It is a one-step
solver – in computing y(tn), it needs
only the solution at the immediately
preceding time point, y(tn-1). In
general, ode45 is the best function to
apply as a first try for most
problems.

ode23 is an implementation of an
explicit Runge-Kutta (2,3) pair of
Bogacki and Shampine. It may be more
efficient than ode45 at crude
tolerances and in the presence of
moderate stiffness. Like ode45, ode23
is a one-step solver.

ode113 is a variable order
Adams-Bashforth-Moulton PECE solver.
It may be more efficient than ode45 at
stringent tolerances and when the ODE
file function is particularly
expensive to evaluate. ode113 is a
multistep solver — it normally needs
the solutions at several preceding
time points to compute the current
solution.

The above algorithms are intended to
solve nonstiff systems. If they appear
to be unduly slow, try using one of
the stiff solvers below.

ode15s is a variable order solver
based on the numerical differentiation
formulas (NDFs). Optionally, it uses
the backward differentiation formulas
(BDFs, also known as Gear's method)
that are usually less efficient. Like
ode113, ode15s is a multistep solver.
Try ode15s when ode45 fails, or is
very inefficient, and you suspect that
the problem is stiff, or when solving
a differential-algebraic problem.

ode23s is based on a modified
Rosenbrock formula of order 2. Because
it is a one-step solver, it may be
more efficient than ode15s at crude
tolerances. It can solve some kinds of
stiff problems for which ode15s is not
effective.

ode23t is an implementation of the
trapezoidal rule using a "free"
interpolant. Use this solver if the
problem is only moderately stiff and
you need a solution without numerical
damping. ode23t can solve DAEs.

ode23tb is an implementation of
TR-BDF2, an implicit Runge-Kutta
formula with a first stage that is a
trapezoidal rule step and a second
stage that is a backward
differentiation formula of order two.
By construction, the same iteration
matrix is used in evaluating both
stages. Like ode23s, this solver may
be more efficient than ode15s at crude
tolerances.

長街聽風 2024-09-10 15:11:22

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

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 

I think this code could work. Try this.

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 
画▽骨i 2024-09-10 15:11:22

代码没问题。只是你必须在 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

我检查了虚拟函数,结果很有希望。

The code is fine. Just you have to add another loop within the for loop. To check the level of consistency.

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

I checked for a dummy function and the results were promising.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文