naeulerb.m

Matlab数值解算法实现代码

function [x,y]=naeulerb(dyfun,xspan,y0,h)
%用途: 隐式Euler法解常微分方程 y'=f(x,y),y(x0)=y0
%格式: [x,y]=naeulerb(dyfun,xspan,y0,h),dyfun为函数f(x,y)，xspan为求解区间[x0,xN]
%      y0为初值y(x0)，h为步长,x返回节点，y返回数值解
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
    y(n+1)=iter(dyfun,x(n+1),y(n),h);
end
x=x';y=y';
function y=iter(dyfun,x,y,h)
y0=y;e=1e-4;K=1e+4;
y=y+h*feval(dyfun,x,y);
y1=y+2*e;k=1;
while abs(y-y1)>e
    y1=y;
    y=y0+h*feval(dyfun,x,y);
    k=k+1;if k>K,error('迭代发散');end
end