myeuler.m

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

function myeuler(f,u0,a,b,h,U)
%EULER法
%调用格式:myeuler(function_handle,u0,a,b,h,function_handle)
%test1: du/dt=u/2;u0=1;h=0.1;a=0;b=1;
%调用:先定义微分方程f=@(t,u) u/2;U=@(t) exp(t/2);myeuler(f,1,0,1,0.1,U)
%test2:典型的伯努利方程：y'=y-2x/y 
%初值为x=0,y=1

%2007-03-08

tic;
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*f(t(i-1),u(i-1));
end
plot(t,u,'b.');
hold on;
line(t,u,'Color','b');
if nargin==6,
	fplot(U,[0,1],'r');
	error=sqrt(sum(abs(U(t)-u).^2))
else
		warning('参数中没有给出原函数！');
end;
toc