teuler.m

来自「求解一阶常微分方程的两个欧拉法」· M 代码 · 共 34 行

function teuler(f,u0,a,b,h,error)
%改进EULER法
%调用格式:teuler(function_handle,u0,a,b,h)
%test: du/dt=u/2;u0=1;h=0.1;a=0;b=1;error=1e-3;
%调用:先定义微分方程f=@(t,u) u/2;teuler(f,1,0,1,0.1,1e-3)

%2007-03-08

tic;
if (nargin<=5),
    error=1e-3;
end;
len=1+(b-a)/h;
u=zeros(1,len);
t=a:h:b;
u(1)=u0;
for i=2:len,
    u(i)=u(i-1)+h*u(i-1)/2;
    temp=u(i)+2*error;
    while abs(u(i)-temp)>error,
        u(i)=u(i-1)+h*(f(t(i-1),u(i-1))+f(t(i),temp))/2;
        temp1=temp;
        temp=u(i);
        u(i)=temp1;
    end
end
plot(t,u,'b+');
hold on;
line(t,u,'Color','b');
fun=@(t) exp(t/2);
fplot(fun,[0,1],'r');
disp('总误差是:');
error=sqrt(sum(abs(fun(t)-u).^2))
toc