计算方法实验报告.docx

计算方法实验报告.docx

《计算方法实验报告.docx》由会员分享，可在线阅读，更多相关《计算方法实验报告.docx（15页珍藏版）》请在冰豆网上搜索。

计算方法实验报告

专业：

电气工程及其自动化

姓名：

学号：

日期：

实验报告

课程名称：

计算方法指导老师：

一、

列主元高斯消去法

源程序（其中方程为,s为主元素绝对值的最小值）：

functionx=gauss_lie（A,b,s）

n=length（b）;

L=A;

b=b';

fork=1:

n-1

[c,i]=max（abs（L（k:

n,k）））;

ik=i+k-1;

ifik~=k

m=L（k,:

）;L（k,:

）=L（ik,:

）;L（ik,:

）=m;

p=b（k）;b（k,:

）=b（ik）;b（ik）=p;

end

ifabs（L（k,k））

disp（'detA=0'）;

return;

end

L（k+1:

n,k）=L（k+1:

n,k）/L（k,k）;

L（k+1:

n,k+1:

n）=L（k+1:

n,k+1:

n）-L（k+1:

n,k）*L（k,k+1:

n）;

b（k+1:

n）=b（k+1:

n）-L（k+1:

n,k）*b（k）;

end

ifabs（L（n,n））

disp（'detA=0'）;

return;

end

b（n）=b（n）/L（n,n）;

fork=n-1:

-1:

1

b（k）=（b（k）-L（k,k+1:

n）*b（k+1:

n））/L（k,k）;

end

x=b;

流程图:

结果显示：

输出方程解

二、

列主元三角分解法

源程序（其中方程为,s为主元素绝对值的最小值）：

functionX=PLU（A,b,s）

n=length（b）;

Lt=A;

b=b';

p=eye（n）;

fork=1:

n-1

[c,i]=max（abs（Lt（k:

n,k）））;

ik=i+k-1;

ifik~=k

m=Lt（k,:

）;Lt（k,:

）=Lt（ik,:

）;Lt（ik,:

）=m;

t=p（k,:

）;p（k,:

）=p（ik,:

）;p（ik,:

）=t;

end

ifabs（Lt（k,k））

disp（'detA=0'）;

return;

end

Lt（k+1:

n,k）=Lt（k+1:

n,k）/Lt（k,k）;

Lt（k+1:

n,k+1:

n）=Lt（k+1:

n,k+1:

n）-Lt（k+1:

n,k）*Lt（k,k+1:

n）;

end

ifabs（Lt（n,n））

disp（'detA=0'）;

return;

end

b=p*b;

b

（1）=b

（1）;

fork=2:

n

b（k）=b（k）-Lt（k,1:

k-1）*b（1:

k-1）;

end

b（n）=b（n）/Lt（n,n）;

fork=n-1:

-1:

1

b（k）=（b（k）-Lt（k,k+1:

n）*b（k+1:

n））/Lt（k,k）;

end

X=b;

P=p

L=eye（n）;

fork=1:

n-1

L（k+1:

n,k）=Lt（k+1:

n,k）;

end

L

fork=1:

n

U（k,k:

n）=Lt（k,k:

n）;

end

U

流程图:

结果显示：

>>A=[0,3,4;1,-1,1;2,1,2];b=[1,2,3];s=0.01;

PLU（A,b,s）

P=

001

100

010

L=

1.000000

01.00000

0.5000-0.50001.0000

U=

212

034

002

ans=

1.1667

-0.3333

0.5000

输出PLU矩阵和方程解

三、最小二乘法拟合多项式

源程序：

functionercheng（x,y,a）

n=length（x）;

fori=1:

a+1

forj=1:

n

c（j,i）=x（j）^（i-1）;

end

end

A=c'*c;

b=c'*y';

ans=A\b;

symsX

Y=0;

fori=1:

a+1

Y=Y+ans（i）*X^（i-1）;

end

Y

k=fliplr（ans'）;

x1=linspace（min（x）,max（x））;

y1=polyval（k,x1）;

plot（x,y,'*',x1,y1）;

s=0;

Max=0;

fori=1:

n

t=0;

forj=1:

a+1

t=t+ans（j）*x（i）^（j-1）;

end

t=abs（y（i）-t）;

ift>Max

Max=t;

end

s=s+t^2;

end

s=sqrt（s）

Max

流程图：

结果显示：

>>x=[2468];

y=[2112840];

ercheng（x,y,1）

Y=

-25/2+131/20*X

s=

3.2711

Max=

2.7000

输出拟合多项式均方误差拟合误差和拟合函数和输入点的图象

>>x=[2468];

y=[2112840];

ercheng（x,y,2）

Y=

-35/4+187/40*X+3/16*X^2

s=

2.9069

Max=

1.9500

输出拟合多项式均方误差拟合误差和拟合函数和输入点的图象

四、龙贝格法

源程序:

（输入函数f,积分起始值,终止值,误差要求,最小二分次数）

functionromberg（f,a,b,e,kmin）

formatlong

t（1,1）=（b-a）/2*（f（a）+f（b））;

k=1;

whilek>0

s=0;

fori=1:

2^（k-1）

s=s+f（a+（2*i-1）*（b-a）/（2^k））;

end

t（k+1,1）=0.5*t（k,1）+（b-a）/（2^k）*s;

ifk<4

fori=1:

k

t（k+1,i+1）=（4^i*t（k+1,i）-t（k,i））/（4^i-1）;

end

ifabs（t（k+1,k+1）-t（k,k））

ifk>=kmin

ans=t（k+1,k+1）

break;

end

end

else

fori=1:

3

t（k+1,i+1）=（4^i*t（k+1,i）-t（k,i））/（4^i-1）;

end

ifabs（t（k+1,4）-t（k,4））

ifk>=kmin

ans=t（k+1,4）

break;

end

end

end

k=k+1;

end

t

流程图：

结果显示：

对最小二分次数不做规定

>>f=inline（'sin（x）/x','x'）

f=

Inlinefunction:

f（x）=sin（x）/x

>>romberg（f,10^（-100）,1,5*10^（-7）,1）

ans=

0.94608307038722

t=

0.92073549240395000

0.939793284806180.9461458822735900

0.944513521665390.946086933951790.946083004063670

0.945690863582700.946083310888470.946083069350920.94608307038722

输出最后结果和龙贝格矩阵

规定最小二分次数为6次

>>romberg（f,10^（-100）,1,5*10^（-7）,6）

ans=

0.94608307036718

t=

0.92073549240395000

0.939793284806180.9461458822735900

0.944513521665390.946086933951790.946083004063670

0.945690863582700.946083310888470.946083069350920.94608307038722

0.945985029934390.946083085384950.946083070351380.94608307036726

0.946058560962770.946083071305560.946083070366940.94608307036718

0.946076943060060.946083070425830.946083070367180.94608307036718