poster
2019-12-10, 20:41
这是我使用后向Euler (http://en.wikipedia.org/wiki/Backward_Euler_method)方法以数值方式求解ODE (http://en.wikipedia.org/wiki/Ordinary_differential_equation)的MATLAB (http://en.wikipedia.org/wiki/MATLAB) / FreeMat (http://en.wikipedia.org/wiki/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向后欧拉solve in y1: y1 - h*f(x1,y1) = y0
您的方法： y1 = y0 +h*f(x0,x0+h*f(x0,y0))
您的方法不落后于欧拉。
您无需在y1求解，只需使用正向Euler方法估算y1 。我不想对您的方法进行分析，但是我相信，即使与正向Euler相比，它的表现也确实很差，因为您在错误的点评估了函数f 。
这也是我能想到的最接近您的方法的方法，也是显式的，它应能提供更好的结果。这是Heun的方法 (http://en.wikipedia.org/wiki/Heun's_method) ：
y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))
更多&回答... (https://stackoverflow.com/questions/2937183)
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向后欧拉solve in y1: y1 - h*f(x1,y1) = y0
您的方法： y1 = y0 +h*f(x0,x0+h*f(x0,y0))
您的方法不落后于欧拉。
您无需在y1求解，只需使用正向Euler方法估算y1 。我不想对您的方法进行分析，但是我相信，即使与正向Euler相比，它的表现也确实很差，因为您在错误的点评估了函数f 。
这也是我能想到的最接近您的方法的方法，也是显式的，它应能提供更好的结果。这是Heun的方法 (http://en.wikipedia.org/wiki/Heun's_method) ：
y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))
更多&回答... (https://stackoverflow.com/questions/2937183)