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

10.6.4 米尔恩公式及其MATLAB程序

10.6.4  米尔恩公式及其MATLAB程序

用米尔恩(Milne)公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

function [k,X,Y,wucha,P]= Milne(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)=Y(k-3)+(4*h/3)*((f(k-2:k))'*[2 -1 2]');

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.4  先利用常用的四阶龙格-库塔公式求初值问题的几个数值解,再分别利用米尔恩公式和改进的亚当斯方法及常用的四阶龙格-库塔公式求解其余的数值解,h=1/4,并计算它们与精确解的误差,在同一图形窗口画出精确解和数值解的图形.

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

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

[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;

y=X-1+exp(-X); 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,Y,'gh',X,Y1,'mp',X,Y2,'ro',X,y,'b-'),

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

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

wchY=abs(y-Y),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,Y,Y1,Y2,wchY,wchY1,wchY2],

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

用米尔恩方法求解常微分方程初值问题(10.5)的数值解的MATLAB主程序

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

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

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

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

Y1(k+1)=Y(k-1)+(h/3)*((f(k-1:k+1))'*[1 4 1]');

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

Y(k)=Y1(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.5  先利用常用的四阶龙格-库塔公式求初值问题的几个数值解,再分别利用米尔恩方法,米尔恩公式和改进的亚当斯方法及常用的四阶龙格-库塔公式求解其余的数值解,h=1/4,并计算它们与精确解的误差,在同一图形窗口画出精确解和数值解的图形.

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

>>x0=0;b=8;y0=0;h=1/4;

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

[k,X,Yf,wucha,P]= Milne1(@funfcn1,x0,b,y0,h)

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

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

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

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(@funfcn1,@fun1,x0,b,C,y0,h)

plot(X,Y,'mh',X,Yf,'b>',X,Y1,'gp',X,Y2,'ro',X,y,'b-'),

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

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

wchY=abs(y-Y),wchYf=abs(y-Yf),

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,Y,Yf,Y1,Y2,wchY,wchYf,wchY1,wchY2],

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

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