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

10.6.3 亚当斯隐式公式及其MATLAB程序

10.6.3  亚当斯隐式公式及其MATLAB程序

 

例10.6.2  先利用常用的四阶龙格-库塔公式求初值问题

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

解  编写并保存名为Adams4y.m的MATLAB程序

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

x1=x+h/2; x2=x+h/2;

x3=x+h;k1=x-y;

y1=y+h*k1/2;

x=x+h;k2=x1-y1;

y2=y+h*k2/2;k3=x2-y2;

y3=y+h*k3;k4=x3-y3;

y=y+h*(k1+2*k2+2*k3+k4)/6;

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

end

X,Y,

for k=3:n

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

Y(k+1)=(1/24.9)*(0.24*k+0.12+(Y(k-2:k))'*[-0.1 0.5 22.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'];

 

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

>> x0=0;b=1;y0=0;h=1/10;

[k,X,Y,wucha,P]= Adams4y (x0,b,y0,h)

y=X-1+exp(-X);

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=x-y,y(0)=0在[0,1]上的数值解','用常用的四阶龙格-库塔公式计算dy/dx=x-y,y(0)=0在[0,1]上的数值解',' dy/dx=x-y,y(0)=0的精确解y=x-1+exp(-x)')

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],

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

 

(二)改进的亚当斯方法及其MATLAB程序

用改进的亚当斯方法求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

function [k,X,Y,wucha,P]=Adams4y1(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);x=x+h;

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

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

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

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

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.3  用改进的亚当斯方法和常用的四阶龙格-库塔公式求例10.6.2中初值问题的数值解及其误差,在同一图形窗口画出精确解和数值解的图形,并与例10.6.2的图形比较.

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

>> x0=0;b=1;

y0=0;h=1/10;

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

y=X-1+exp(-X);

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,'ro'),

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

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

wchY=abs(y-Y),

wchY1=abs(y-Y1),

m=zeros(1,k),

for n=1:k,

m(1,n)=n,

end, [

m',X,y,Y,Y1,wchY,wchY1],

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

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