向前欧拉公式计算matlab程序函数

向前欧拉公式计算matlab程序函数

function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol)

%初始化.

pow=1/3;

if nargin < 5 | isempty(tol),

tol = 1.e-6;

end;

x=x0;

h=0.0078125*(b-x);

y=y0(:);

p=128;

H=zeros(p,1);

X=zeros(p,1);

Y=zeros(p,length(y));

k=1;

X(k)=x;

Y(k,:)=y';

% 绘图.

clc,x,h,y

% end

% 主循环.

while (xx)

if x+h>b

h=b-x;

end

%计算斜率.

fxy=feval(funfcn,x,y);

fxy=fxy(:);

%计算误差，设定可接受误差.

delta=norm(h*fxy,'inf');

wucha=tol*max(norm(y,'inf'),1.0);

% 当误差可接受时重写解.

if delta<=wucha

x=x+h; y=y+h*fxy; k=k+1;

if k>length(X)

X=[X;zeros(p,1)];

Y=[Y;zeros(p,length(y))];

H=[H;zeros(p,1)];

end

H(k)=h;

X(k)=x;

Y(k,:)=y';

plot(X,Y,'rp'),

grid

xlabel('自变量X'),

ylabel('因变量Y')

end

%更新步长.

if delta~=0.0

h=min(h*8,0.9*h*(wucha/delta)^pow);

end

end

if (x

disp('Singularity likely.'), x

end

H=H(1:k);

X=X(1:k);

Y=Y(1:k,:);

n=1:k;

P=[n',H,X,Y]

end

function y=funfcn(x,y)

y(1)=-3*y+8*x-7;

end

clc,clear ,close all

subplot(2,1,1)

x0=0;

y0=1;

b=2;

n=10;

h=2/10;

[k,X,Y,P,REn]=QEuler(@funfcn,x0,y0,b,h)

hold on

S1=1+1/6*(6+12*X+30*exp(2*X)).^(1/2)

plot(X,S1,'b-')

title('用向前欧拉公式计算dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解') legend('n=10数值解','精确解')

hold off

jdwucY=S1-Y;

jwY=S1-Y;

xwY=jwY./Y;

k1=1:n;

k=[0,k1];

% P1=[k',X,Y,S1,jwY,xwY]

subplot(2,1,2)

n1=100;

h1=2/100;

[k,X1,Y1,P1,Ren1]=QEuler(@funfcn,x0,y0,b,h1) hold on

S2=1+1/6*(6+12*X1+30*exp(2*X1)).^(1/2) plot(X1,S2,'b-')

legend('n=100数值解','精确解')

hold off

jwY1=S2-Y1;

xwY1=jwY1./Y1;

k1=1:n1;k=[0,k1];

% P2=[k',X1,Y1,S2,jwY1,xwY1]