数值分析报告编程及运行结果高斯顺序消元法.docx

数值分析报告编程及运行结果高斯顺序消元法.docx

《数值分析报告编程及运行结果高斯顺序消元法.docx》由会员分享，可在线阅读，更多相关《数值分析报告编程及运行结果高斯顺序消元法.docx（24页珍藏版）》请在冰豆网上搜索。

数值分析报告编程及运行结果高斯顺序消元法

高斯消元法

1.程序:

clear

formatrat

A=input（'输入增广矩阵A='）

[m,n]=size（A）;

fori=1:

（m-1）

numb=int2str（i）;

disp（['第',numb,'次消元后的增广矩阵']）

forj=（i+1）:

m

A（j,:

）=A（j,:

）-A（i,:

）*A（j,i）/A（i,i）;

end

A

end

%回代过程

disp（'回代求解'）

x（m）=A（m,n）/A（m,m）;

fori=（m-1）:

-1:

1

x（i）=（A（i,n）-A（i,i+1:

m）*x（i+1:

m）'）/A（i,i）;

end

x

2.运行结果：

高斯选列主元消元法

1.程序:

clear

formatrat

A=input（'输入增广矩阵A='）

[m,n]=size（A）;

fori=1:

（m-1）

numb=int2str（i）;

disp（['第',numb,'次选列主元后的增广矩阵']）

temp=max（abs（A（i:

m,i）））;

[a,b]=find（abs（A（i:

m,i））==temp）;

tempo=A（a

（1）+i-1,:

）;

A（a

（1）+i-1,:

）=A（i,:

）;

A（i,:

）=tempo

disp（['第',numb,'次消元后的增广矩阵']）

forj=（i+1）:

m

A（j,:

）=A（j,:

）-A（i,:

）*A（j,i）/A（i,i）;

end

A

end

%回代过程

disp（'回代求解'）

x（m）=A（m,n）/A（m,m）;

fori=（m-1）:

-1:

1

x（i）=（A（i,n）-A（i,i+1:

m）*x（i+1:

m）'）/A（i,i）;

end

x

2.运行结果：

追赶法

1.程序：

function[x,L,U]=zhuiganfa（a,b,c,f）

a=input（'输入矩阵-1对角元素a='）;

b=input（'输入矩阵对角元素b='）;

c=input（'输入矩阵+1对角元素c='）;

f=input（'输入增广矩阵最后一列元素f='）;

n=length（b）;

%对A进行分解

u

（1）=b

（1）;

fori=2:

n

if（u（i-1）~=0）

l（i-1）=a（i-1）/u（i-1）;

u（i）=b（i）-l（i-1）*c（i-1）;

else

break;

end

end

L=eye（n）+diag（l,-1）;

U=diag（u）+diag（c,1）;

x=zeros（n,1）;

y=x;

%求解Ly=b

y

（1）=f

（1）;

fori=2:

n

y（i）=f（i）-l（i-1）*y（i-1）;

end

%求解Ux=y

if（u（n）~=0）

x（n）=y（n）/u（n）;

end

fori=n-1:

-1:

1

x（i）=（y（i）-c（i）*x（i+1））/u（i）;

end

2.运行结果：

高斯-塞德尔迭代格式

1.程序：

functionx=Gauss_Seidel（a,b）

a=input（'输入系数矩阵a='）

b=input（'输入增广矩阵最后一列b='）;

e=0.5e-7;

n=length（b）;

N=50;

x=zeros（n,1）;

t=zeros（n,1）;

fork=1:

N

sum=0;

E=0;

t（1:

n）=x（1:

n）;

fori=1:

n

x（i）=（b（i）-a（i,1:

（i-1））*x（1:

（i-1））-a（i,（i+1）:

n）*t（（i+1）:

n））/a（i,i）;

end

ifnorm（x-t）

k

break;

end

end

2.运行结果：

雅戈比迭代格式

1.程序：

functionx=Jocabi（a,b）

a=input（'输入系数矩阵a='）;

b=input（'输入增广矩阵最后一列b='）;

e=0.5e-7;

n=length（b）;

N=100;

x=zeros（n,1）;

y=zeros（n,1）;

fork=1:

N

sum=0;

fori=1:

n

y（i）=（b（i）-a（i,1:

n）*x（1:

n）+a（i,i）*x（i））/a（i,i）;

end

fori=1:

n

sum=sum+（y（i）-x（i））^2;

end

ifsqrt（sum）

k

break;

else

fori=1:

n

x（i）=y（i）;

end

end

end

ifk==Nwarning（'未能找到近似解'）;

end

2.运行结果：

逐次超松弛法（SOR）

1.程序：

function[n,x]=sor22（A,b,X,nm,w,ww）

%用超松弛迭代法求解方程组Ax=b

%输入：

A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子

%输出：

x为求得的方程组的解构成的列向量，n为迭代次数

A=input（'输入系数矩阵A='）;

b=input（'输入方程组右端的列向量b='）;

X=input（'输入迭代初值构成的列向量X='）;

nm=input（'输入最大迭代次数nm='）;

w=input（'输入误差精度w='）;

ww=input（'输入松弛因子ww='）;

n=1;

m=length（A）;

D=diag（diag（A））;%令A=D-L-U,计算矩阵D

L=tril（-A）+D;%令A=D-L-U,计算矩阵L

U=triu（-A）+D;%令A=D-L-U,计算矩阵U

M=inv（D-ww*L）*（（1-ww）*D+ww*U）;%计算迭代矩阵

g=ww*inv（D-ww*L）*b;%计算迭代格式中的常数项

%下面是迭代过程

whilen<=nm

x=M*X+g;%用迭代格式进行迭代

ifnorm（x-X,'inf'）

disp（'迭代次数为'）;n

disp（'方程组的解为'）;x

return;

%上面：

达到精度要求就结束程序，输出迭代次数和方程组的解

end

X=x;n=n+1;

end

%下面：

如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）

disp（'在最大迭代次数内不收敛!

'）;

disp（'最大迭代次数后的结果为'）;

2.运行结果：

二分法求解方程的根

1.程序：

%其中a，b表示查找根存在的范围，M表示要求解函数的值

%f（x）表示要求解根的方程

%eps表示所允许的误差大小

functiony=er_fen_fa（a,b,M）

k=0;

eps=0.05

whileb-a>eps

x=（a+b）/2;

%检查是否大于值

if（（x^3）-3*x-1）>M

b=x

else

a=x

end

k=k+1

end

2.运行结果：

Newton迭代法（切线法）

1.程序：

functionx=nanewton（fname,dfname,x0,e,N）

%newton迭代法解方程组

%fname和dfname分别表示F（x）及其导函数的M函数句柄或内嵌函数，x0为迭代初值，e为精度要求

x=x0;x0=x+2*e;k=0;

ifnargin<5,N=500;end

ifnargin<4e=1e-4;end

whileabs（x0-x）>e&k

k=k+1;

x0=x;x=x0-feval（fname,x0）/feval（dfname,x0）;

disp（x）

end

ifk==N,warning（'已达迭代次数上限'）;

end

2.运行结果：

割线方式迭代法

1.程序：

functionx=ge_xian_fa（fname,dfname,x0,x1,e,N）

%割线方式迭代法解方程组

%fname和dfname分别表示F（x）及其导函数的M函数句柄或内嵌函数，x0,x1分别为迭代初值，e为精度要求

k=0;a=x1;b=x0;

ifnargin<5,N=500;end

ifnargin<4e=1e-4;end

whileabs（a-b）>e&k

k=k+1;

x=x1-（（x1-x0）/（feval（fname,x1）-feval（fname,x0）））*feval（fname,x1）;

iffeval（fname,x）*feval（fname,x0）>0,

x0=x;b=x0;

else

x1=x;a=x1;

end

x=x1-（（x1-x0）/（feval（fname,x1）-feval（fname,x0）））*feval（fname,x1）;

numb=int2str（k）;

disp（['第',numb,'次计算后x=']）

fprintf（'%f\n\n',x）;

end

ifk==N,warning（'已达迭代次数上限'）;

end

2.运行结果：

Newton插值

1.程序：

%保存文件名为New_Int.m

%Newton基本插值公式

%x为向量，全部的插值节点

%y为向量，差值节点处的函数值

%xi为标量，是自变量

%yi为xi出的函数估计值

functionyi=newton_chazhi（x,y,xi）

n=length（x）;

m=length（y）;

ifn~=m

error（'ThelengthsofXangYmustbeequal!

'）;

return;

end

%计算均差表Y

Y=zeros（n）;

Y（:

1）=y';

fork=1:

n-1

fori=1:

n-k

ifabs（x（i+k）-x（i））

error（'theDATAiserror!

'）;

return;

end

Y（i,k+1）=（Y（i+1,k）-Y（i,k））/（x（i+k）-x（i））;

end

end

%计算牛顿插值公式

yi=0;

fori=1:

n

z=1;

fork=1:

i-1

z=z*（xi-x（k））;

end

yi=yi+Y（1,i）*z;

end

2.运行结果：

Lagrange插值

1.程序：

functiony0=Language（x,y,x0）

symstl;

iflength（x）==length（y）

n=length（x）;

else

disp（'x和y的维数不相等！

'）;

return;%检错

end

h=sym（0）;

fori=1:

n

l=sym（y（i））;

forj=1:

i-1

l=l*（t-x（j））/（x（i）-x（j））;

end;

forj=i+1:

n

l=l*（t-x（j））/（x（i）-x（j））;

end;

h=h+l;

end

simplify（h）;

ifnargin==3

y0=subs（h,'t',x0）;%计算插值点的函数值

else

y0=collect（h）;

y0=vpa（y0,6）;%将插值多项式的系数化成6位精度的小数

end

2.运行结果：

最小二乘法

1.程序：

functionp=nafit（x,y,m）

%多项式拟合

%x,y为已知数据点向量,分别表示横,纵坐标,m为拟合多项式的次数,结果返回m次拟合多项式系数,从高次到低次存放在向量p中.

A=zeros（m+1,m+1）;

fori=0:

m

forj=0:

m

A（i+1,j+1）=sum（x.^（i+j））;

end

b（i+1）=sum（x.^i.*y）;

end

a=A\b';

p=fliplr（a'）;

t=0:

0.1:

1.6;

S3=polyval（p,t）;

plot（x,y,'pk'）;

holdon

plot（t,S3,'r'）;

2.运行结果：