f=inline('cos(x)+sin(y)','x','y'); %微分方程的右边项
dx=0.05; %x方向步长
xleft=pi/2; %区域的左边界
xright=3*pi/2; %区域的右边界
xx=xleft:dx:xright; %一系列离散的点
n=length(xx); %点的个数
y0=0;
%%(1)欧拉法
Euler=y0;
for i=2:n
Euler(i)=Euler(i-1)+dx*f(xx(i-1),Euler(i-1));
end
%%(2)改进欧拉法
MEuler=y0;
for i=2:n
MEuler(i)=MEuler(i-1)+dx/2*(f(xx(i-1),MEuler(i-1))+f(xx(i),MEuler(i-1)+dx*f(xx(i-1),MEuler(i-1))));
end
%%(3)龙格库塔法
RK=y0;
for i=2:n
k1=f(xx(i-1),RK(i-1));
k2=f(xx(i-1)+dx/2,RK(i-1)+k1*dx/2);
k3=f(xx(i-1)+dx/2,RK(i-1)+k2*dx/2);
k4=f(xx(i-1)+dx,RK(i-1)+k3*dx);
RK(i)=RK(i-1)+dx*(k1+2*k2+2*k3+k4)/6;
end
%%Euler和MEuler,RK法结果比较
plot(xx,Euler,xx,MEuler,xx,RK)
hold on
legend('Euler','MEuler','Runge-Kutta')