有限差分法Eluer算法（求解常微分方程）

在MATLAB 里面写一个程序：要求用 隐式欧拉法(backward euler) 去解决常微分方程。
下面是两个例题，
1. x' = - 2x ; 还给出准确值是： x= e^(-2t) 要求： 求出t=0， 这是一个初始值，然后算出在区间[0,5] 的值。
还给出 步长(step size) h=0.1
2. x' = xsint - 2x^2; 给出初始条件 x(0) = 0
要求： 求出在区间[0,1] 的值，比较 h=0.01 和 h=0.001 在 t = 1 这里。

1.新建一个m文件，编写隐式Euler法的程序：
function [x,y]=Implicit_Euler(odefun,xspan,y0,h,varargin)
% 隐式Euler公式求解常微分方程
% 输入参数：
% ---odefun：微分方程的函数描述
% ---xspan：求解区间[x0,xn]
% ---y0：初始条件
% ---h：迭代步长
% ---p1,p2,…：odefun函数的附加参数
% 输出参数：
% ---x：返回的节点，即x=xspan(1):h:xspan(2)
% ---y：微分方程的数值解
x=xspan(1):h:xspan(2);
y(1)=y0;
for k=1:length(x)-1
z0=y(k)+h*feval(odefun,x(k),y(k),varargin{:});
z1=inf;
while abs(z1-z0)>1e-4
z1=y(k)+h*feval(odefun,x(k+1),z0,varargin{:});
z0=z1;
end
y(k+1)=z1;
end
x=x;y=y;

2.在命令窗口直接调用上面的程序，求解常微分问题：
（1）
f=@(t,x)-2*x;
[t,x]=Implicit_Euler(f,[0 5],1,0.1);
t1=0:0.1:5;
x1=exp(-2*t);
plot(t1,x1)
hold on
plot(t,x,'k.-','markersize',16)
legend('解析解','隐式Euler求解结果')
xlabel('t');ylabel('x');

（2）此题给出的初值好像有问题吧，x0=0的话，求解的结果都是为0，所以我改用x0=1求解试了一下：
f=@(t,x)x*sin(t)-2*x^2;
[t,x]=Implicit_Euler(f,[0 1],1,0.01);
[t1,x1]=Implicit_Euler(f,[0 1],1,0.001);
e=x1(end)-x(end)
>>e =
-0.0029

plot(t,x,'r:')
hold on
plot(t1,x1,'g--')
xlabel('t');ylabel('x')
legend('积分步长为0.01','积分步长为0.001')

利用Matlab实现用：
向前Euler方法
向后Euler方法
改进的Euler方法
梯形公式

求方程
y'=y-2x/y,   0<x<1
y(0）=1,
的数值解{yn}
在命令窗口的程序：
f=@(x,y)y-2*x/y;
[x,y]=euler_xh(f,[0 5],1,0.1);
x1=0:0.1:5;
y1=0:0.1:5;
plot(x1,y1,'r.-','markersize',16)
hold on
plot(x,y,'k.-','markersize',16)
legend('解析解','隐式Euler求解结果')
xlabel('t');ylabel('x');

楼主您好，我给您个前向欧拉法和改进欧拉法的
%  forward euler method
function  [ x, y ] = eulerf ( f, y0, a, b, n )
y ( 1 ) = y0 ;
h = ( b - a ) / n;
x = a : h : b;
for i = 1 : n
y ( i +1 ) = y ( i ) + h * feval ( f, x ( i ), y ( i ) ) ;
end;

% Improved euler method
function  [ x, y ] = euleri ( f, y0, a, b, n )
y ( 1 ) = y0;
h = ( b - a ) / n;
x = a : h : b;
for i = 1 : n
yt = y ( i ) + h * feval ( f, x ( i ), y ( i ) );
yi = y ( i ) + h * feval ( f, x ( i + 1 ), yt );
y ( i + 1 ) = 1 / 2 * ( yt + yi );
end;

至于后向欧拉法和梯形公式法，我建议针对f(x,y)的实际函数进行编译，在每个循环步骤中都需要求解方程，
解出y(i+1)的值，但是有时候y(i+1)需要用数值解法来求解，即使可以用solve函数，也会比较慢。。。
如果编辑这两种方法的程序，我再给您贴个后向欧拉法的，需要中间预测y值：
% backward euler method
function  [ x, y ] = eulerb ( f, y0, a, b, n )
y ( 1 ) = y0;
h = ( b - a ) / n;
x = a : h : b;
for i = 1 : n
yt = y ( i ) + h * feval ( f, x ( i ), y ( i ) );
y ( i +1 ) = y ( i ) + h * feval ( f, x ( i ), yt );
end;

下面给您贴两个循环体内迭代求解y(i+1)值的后向欧拉法和梯形公式的程序：

%  Backward euler method
function  [ x, y ] = trapzd ( f, y0, a, b, n )
y ( 1 ) = y0;
h = ( b - a ) / n;
x = a : h : b;
for i = 1 : n
yt = y ( i ) + h * feval ( f, x ( i ), y ( i ) );
done = 0;
while  ~done
y ( i + 1 ) = y ( i ) + h * feval ( f, x ( i ), yt );
done = ( abs ( y ( i + 1 ) - yt ) < 1e-6 );
yt = y ( i + 1 );
end;
end;

%   trapezoidal method
function  [ x, y ] = trapzd ( f, y0, a, b, n )
y ( 1 ) = y0;
h = ( b - a ) / n;
x = a : h : b;
for i = 1 : n
yt = y ( i ) + h * feval ( f, x ( i ), y ( i ) );
done = 0;
while  ~done
y ( i + 1 ) = y ( i ) + 0.5 * h * ( yt + feval ( f, x ( i ), yt ) );
done = ( abs ( y ( i + 1 ) - yt ) < 1e-6 );
yt = y ( i + 1 );
end;
end;

额 已经搞定了 不过还是谢谢 顺便附上向后欧拉和梯形公式
function [x,y]=euler_xh(f,r,y0,h) %f为求解方程,r为x范围,y初值y0,h为步长
x=r(1):h:r(2);
y0=1;
y(1)=y0;
for n=1:length(x)-1
y(n+1)=y(n)+h*feval(f,x(n),y(n));
y(n+1)=y(n)+h*feval(f,x(n+1),y(n+1));
end
x=x';
y=y';

function [x,y]=Traprl(f,r,y0,h) %f为求解方程,r为x范围,y初值y0,h为步长
x=r(1):h:r(2);
y0=1;
y(1)=y0;
for n=1:length(x)-1
y(n+1)=y(n)+0.01*feval(f,x(n),y(n));
y(n+1)=y(n)+0.005*(feval(f,x(n),y(n))+feval(f,x(n+1),y(n+1)));
end
y=y'