MATLAB 代码帮助。后向欧拉法
这是 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 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(6)
您的方法是一种新型方法。它既不是向后的,也不是向前的欧拉。 :-)
前向欧拉:
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 estimatey1
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 functionf
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)))
我能发现的唯一问题是这一行:
应该是:
避免负步数。如果您向负 x 方向移动,请确保为函数指定负步长。
您的答案可能会有所偏差,因为您对答案的近似程度有多粗略。为了获得半准确的结果,deltaX 必须非常非常小,并且步长必须非常非常小。
附言。这不是“向后欧拉方法”,它只是常规的旧欧拉方法。
如果这是家庭作业,请标记它。
The only issue I can spot is that the line:
Should be:
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.
看看数值食谱,特别是第16章,常微分方程的积分。众所周知,欧拉方法存在问题:
稳定除非您知道您的教科书使用的是欧拉方法,否则我不希望结果匹配。即使是这样,您可能也必须使用相同的步长才能获得相同的结果。
Have a look at numerical recipes, specifically chapter 16, integration of ordinary differential equations. Euler's method is known to have problems:
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.
除非您确实想通过自己编写的欧拉方法求解 ODE,否则您应该看看 内置 ODE 求解器。
旁注:您不需要在循环内创建
x(i)
,如下所示:x(i+1) = x(i)+h;
。相反,您可以简单地编写x = xinit:h:xfinal;
。另外,您可能需要编写n = round(xfinal-xinit)/h);
以避免警告。以下是 MATLAB 实现的求解器。
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 writex = xinit:h:xfinal;
. Also, you may want to writen = round(xfinal-xinit)/h);
to avoid warnings.Here are the solvers implemented by MATLAB.
我认为这段代码可以工作。试试这个。
I think this code could work. Try this.
代码没问题。只是你必须在 for 循环中添加另一个循环。检查一致性水平。
我检查了虚拟函数,结果很有希望。
The code is fine. Just you have to add another loop within the for loop. To check the level of consistency.
I checked for a dummy function and the results were promising.