用自适应梯形公式计算dy

用自适应梯形公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解

%%

clc,clear,close all

x0=0;

y0=1;

b=2;

tol=1.e-1;

[Ht,X,Yt,k,h,Pt]= odtixing2(@funfcn,x0,b,y0,tol),

hold on

S1 = 8/3*X-29/9+38/9*exp(-3*X),

plot(X,S1,'b-'),

hold off

hold on,

[Hq,X,Yq,k,h,Pq]=QEuler(@funfcn,x0,b,y0,tol),

hold off

grid

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

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

%初始化.

pow=1/3;

if nargin < 5 | isempty(tol),

tol = 1.e-6;

end;

% if nargin < 6 | isempty(trace),

% trace = 0;

% end;

x=x0;

h=0.0078125*(b-x);

y=y0(:);

p=128;

n=fix((b-x0)/h);

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

% 主循环.

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; y1=y+h*fxy; fxy1=feval(funfcn,x,y1); fxy=fxy(:);

y2=y+h*fxy1; y=(y1+y2)/2; 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,'go'), grid

xlabel('自变量X'), ylabel('因变量Y')

title('用自适应梯形公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解') 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]