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

10.6.2 亚当斯(Adams)显式公式及其MATLAB程序


10.6.2  亚当斯(Adams)显式公式及其MATLAB程序

用四阶亚当斯显式公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

function [k,X,Y,wucha,P]= Adams4x(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

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; x=x+h;

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;

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)=Y(k)+(h/24)*((f(k-3:k))'*[-9 37 -59 55]');

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)); k=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.1  先利用常用的四阶龙格-库塔公式求初值问题

的几个点的数值解,再利用四阶亚当斯显式公式求解常微分方程初值问题,h=1/15,并计算它与精确解的误差,在同一图形窗口画出精确解和数值解的图形.

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

>> y=dsolve('Dy=1-(2*x*y)/(1+x^2)','x')

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

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

y=(X+1/3*X.^3)./(1+X.^2); b31=0; b41=0;b42=0;b43=1;

c1=1/6;c2=2/6;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,Y1,fxy,wch,wucha,P]=RK4(@funfcn,@fun,x0,b,C,y0,h)

plot(X,Y,'gh',X,Y1,'mp',X,y,'bo'),

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

legend('用四阶亚当斯显式公式计算dy/dx=1-(2xy)/(1+x^2),y(0)=0在[0,2]上的数值解','用常用的四阶龙格-库塔公式计算dy/dx=1-(2xy)/(1+x^2),y(0)=0在[0,2]上的数值解','dy/dx=1-(2xy)/(1+x^2)的精确解y=(x+1/3x^3)/(1+x^2)')

wchY=abs(y-Y), wchY1=abs(y-Y1), m=zeros(1,k),

for n=1:k,m(1,n)=n-1,end, [m',X,y,Y,Y1,wchY,wchY1],

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

 



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