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

10.6.6 预测-校正系统及其MATLAB程序

10.6.6  预测-校正系统及其MATLAB程序

(一) 单环节的亚当斯预测-校正公式及其MATLAB程序

用单环节的亚当斯预测-校正公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

function [k,X,Y,wucha,P]=dAdamsyx(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=feval(funfcn,X(1:4),

Y(1:4));f=f',

for k=4:n

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

P=Y(k)+(h/24)*((f(k-3:k))*[-9 37 -59 55]');

f=[f(2) f(3) f(4) feval(funfcn,X(k+1),P)],

Y(k+1)=Y(k)+(h/24)*(f*[1 -5 19 9]');

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

end

for k=1:n

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

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.7  先利用常用的四阶龙格-库塔公式求例10.6.1的初值问题算出几个点的数值解,再分别利用单环节的亚当斯预测-校正公式、四阶亚当斯显式公式(10.50)、四阶亚当斯隐式公式(10.52)求解常微分方程初值问题,h=1/5,并计算它与精确解y=(x+1/3x3)/(1+x2)的误差,在同一图形窗口画出精确解和数值解的图形.

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

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

subplot(3,1,1)

[k,X,Y1,wucha1,P1]=dAdamsyx(@funfcn,x0,b,y0,h)

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

plot(X,Y1,'mh',X,y,'bo'),

grid

legend('用单环节的亚当斯预测-校正公式计算的数值解','精确解y=(x+1/3x^3)/(1+x^2)')

wch1=abs(y-Y1),

[P1,y,wch1]

title('dy/dx=1-(2xy)/(1+x^2),y(0)=0在[0,2]上的数值解和精确解的图形')

 subplot(3,1,2)

 [k,X,Y2,wucha2,P2]=Adams4x(@funfcn,x0,b,y0,h)

y=(X+1/3*X.^3)./(1+X.^2);plot(X,Y2,'mh',X,y,'bo'),

grid

legend('用四阶亚当斯显式公式计算的数值解','精确解y=(x+1/3x^3)/(1+x^2)')

wch2=abs(y-Y2),

[P2,y,wch2]

subplot(3,1,3)

[k,X,Y,wucha,P]=Adams4y1(@funfcn,x0,b,y0,h),y=(X+1/3*X.^3)./(1+X.^2);

plot(X,Y,'mh',X,y,'bo'),

grid,

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

legend('用四阶亚当斯隐式公式计算的数值解','精确解y=(x+1/3x^3)/(1+x^2)')

wch=abs(y-Y),[P,y,wch],A=[X,y,Y1,Y2,Y,wch1,wch2,wch]

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

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