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

10.6.5 汉明公式及其MATLAB程序

10.6.5  汉明公式及其MATLAB程序

汉明公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

function [k,X,Y,wucha,P]= Hamming(funfcn,x0,b,y0,h)

x=x0;y=y0;p=128; n=fix((b-x0)/h);

if n<5,return,end;

X=zeros(p,1);

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

f=zeros(p,1);

k=1; X(k)=x; Y(k,:)=y';

for k=2:4

x=x+h; c1=1/6;c2=2/6;

c3=2/6;c4=1/6;a2=1/2; a3=1/2;

a4=1;b21=1/2;b31=0;

b32=1/2; b41=0;b42=0;

b43=1;x1=x+a2*h;

x2=x+a3*h; x3=x+a4*h;

k1=feval(funfcn,x,y);

y1=y+b21*h*k1;

k2=feval(funfcn,x1,y1);

y2=y+b31*h*k1+b32*h*k2;

k3=feval(funfcn,x2,y2);

y3=y+b41*h*k1+b42*h*k2+b43*h*k3;

k4=feval(funfcn,x3,y3);

y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4);

X(k)=x;Y(k,:)=y;k=k+1;

end

X;Y; f(1:4)=feval(funfcn,X(1:4),Y(1:4));

for k=4:n

f(k)=feval(funfcn,X(k),Y(k));

X(k+1)=X(1)+h*k;

Y(k+1)=(1/8)*(9*Y(k)-Y(k-2))+(3*h/8)*((f(k-2:k))'*[-1 2 1]');

f(k+1)=feval(funfcn,X(k+1),Y(k+1));

f(k)=f(k+1); k=k+1;

end

for k=2:n+1

wucha(k)=norm(Y(k)-Y(k-1));

end

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

wucha=wucha(1:n,:);P=[n',X,Y,wucha'];

 

例10.6.6  先利用常用四阶龙格-库塔公式求初值问题  的几个点的数值解,再利用汉明公式、米尔恩公式和改进的四阶亚当斯隐式公式及常用的四阶龙格-库塔公式求解其余的数值解求解常微分方程初值问题,h=1/2,并计算它与精确解y=(x+1/3x3)/(1+x2)的误差,在同一图形窗口画出精确解和数值解的图形.

解  在MATLAB工作窗口输入下面的程序

>> x0=0;b=20;y0=0;h=1/2;

[k,X,Yh,wucha,P]= Hamming(@funfcn,x0,b,y0,h)

y=(X+1/3*X.^3)./(1+X.^2);

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

[k,X,Y1,wucha,P]=Adams4y1(@funfcn,x0,b,y0,h),

c1=1/6;c2=2/6;

b31=0; b41=0;b42=0;b43=1;

c3=2/6;c4=1/6;a2=1/2; a3=1/2;

a4=1;b21=1/2; b32=1/2;

C=[c1,c2,c3, c4,a2, a3, a4,b21,b31,b32,b41,b42,b43];

[k,X,Y2,fxy,wch,wucha,P]=RK4(@funfcn,@fun,x0,b,C,y0,h)

plot(X,Yh,'bh',X,Y,'m*',X,Y1,'gp',X,Y2,'ro',X,y,'k-'),

grid

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

legend('用汉明公式计算y=1-2xy/(1+x^2),y(0)=0在[0,20]上的数值解','用米尔恩公式计算y=1-2xy/(1+x^2),y(0)=0在[0,20]上的数值解','用改进的四阶亚当斯方法计算y=1-2xy/(1+x^2),y(0)=0在[0,20]上的数值解','用常用的四阶龙格-库塔公式计算y=1-2xy/(1+x^2),y(0)=0在[0,20]上的数值解',' y=1-2xy/(1+x^2)的精确解y=(x+1/3x^3)/(1+x^2)')

wchY=abs(y-Y),

wchYh=abs(y-Yh),wchY1=abs(y-Y1),

wchY2=abs(y-Y2), m=zeros(1,k),

for n=1:k,

m(1,n)=n-1,

end,

[m',X,y,Yh,Y,Y1,Y2,wchYh,wchY,wchY1,wchY2],

运行后屏幕显示图和计算结果(略).

 

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