matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

10.4.2 梯形公式的两种MATLAB程序

10.4.2  梯形公式的两种MATLAB程序

 

(一)  梯形公式的MATLAB程序

用梯形公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序1

function [X,Y,n,P]= odtixing1(funfcn,x0,b,y0,h,tol)

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

X=zeros(n+1,1);

Y=zeros(n+1,1);

k=1; X(k)=x0;

Y(k,:)=y0; Y1(k,:)=y0;

% 绘图.

clc,x0,h,y0

% 产生初值.

for i=2:n+1

X(i)=x0+h;

fx0y0=feval(funfcn,x0,y0);

Y(i,:)=y0+h*fx0y0;

fxiyi=feval(funfcn,X(i),Y(i,:));

Y1(i,:)=y0+h*(fxiyi+fx0y0)/2;

% 主循环.

Wu=abs(Y1(i,:)-Y(i,:));

while Wu>tol

      p=Y1(i,:),

fxip=feval(funfcn,X(i),p);

Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:), Y(i,:)=p1;

    end

x0=x0+h; y0=Y1(i,:);

Y(i,:)=y0; plot(X,Y,'ro')

grid on

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

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

end

  X=X(1:n+1); Y=Y(1:n+1,:);

n=1:n+1; P=[n',X,Y]

 

例10.4.1  用梯形公式求解区间上的初值问题,取步长,精度为,并与精确解作比较,在同一个坐标系中画出图形.

 

解  输入程序

>> x0=0;y0=1; b=2; tol=0.1; h=0.05;

[X,Yt,n,Pt]=odtixing1(@funfcn,x0,b,y0,h,tol)

hold on

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

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

legend('h=0.05,用梯形公式计算dy/dx=8x-3y-7,y(0)=1在[0,2]上的数值解','dy/dx=8x-3y-7,y(0)=1在[0,2]上的精确解')

juwYt=S1-Yt; xiwYt=juwYt./Yt; Lt=[Pt,S1,juwYt,xiwYt]

运行后屏幕显示取精度为,分别用梯形公式和向前欧拉公式求解此初值问题在区间上的自变量X处数值解Yi(i=t,q)和精确解S1,步长H, Yi的相对误差xiwYi和绝对误差juwYi (略) 及其数值解和精确解的图形.

 

 (二) 自适应梯形公式的MATLAB程序

用自适应梯形公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

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 (x<b)&(x+h>x)

  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<b)

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]

 

例10.4.2  分别用自适应梯形公式和向前欧拉公式分别求解区间上的初值问题,取精度为,并与精确解作比较,在同一个坐标系中作出图形.

解  输入程序

>> 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('用自适应梯形公式计算dy/dx=8x-3y-7,y(0)=1在[0,2]上的数值解')

运行后屏幕显示取精度为,分别用梯形公式和向前欧拉公式求解此初值问题在区间上的自变量X处数值解Yi(i=t,q)和精确解S1及其图形,步长H等 (略).

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3100402.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1