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

高斯消元法1.程序:

clear

format rat

A=input('输入增广矩阵A=')

[m,n]=size(A);

for i=1:(m-1)

numb=int2str(i);

disp(['第',numb,'次消元后的增广矩阵']) for j=(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);

for i=(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

format rat

A=input('输入增广矩阵A=')

[m,n]=size(A);

for i=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,'次消元后的增广矩阵'])

for j=(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);

for i=(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);

for i=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);

for i=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

for i=n-1:-1:1

x(i)=(y(i)-c(i)*x(i+1))/u(i); end

2.运行结果：

高斯-塞德尔迭代格式

1.程序：

function x=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);

for k=1:N

sum=0;

E=0;

t(1:n)=x(1:n);

for i=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

if norm(x-t)

k

break;

end

end

2.运行结果：

雅戈比迭代格式1.程序：

function x=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);

for k=1:N

sum=0;

for i=1:n

y(i)=(b(i)-a(i,1:n)*x(1:n)+a(i,i)*x(i))/a(i,i);

end

for i=1:n

sum=sum+(y(i)-x(i))^2;

end

if sqrt(sum)

k

break;

else

for i=1:n

x(i)=y(i);

end

end

end

if k==N warning('未能找到近似解'); 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; %计算迭代格式中的常数项

%下面是迭代过程

while n<=nm

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

if norm(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表示所允许的误差大小

function y=er_fen_fa(a,b,M)

k=0;

eps=0.05

while b-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.程序：

function x=nanewton(fname,dfname,x0,e,N)

%newton迭代法解方程组

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

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

if nargin<5,N=500;end

if nargin<4 e=1e-4;end

while abs(x0-x)>e&k

k=k+1;

x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);

disp(x)

end

if k==N,warning('已达迭代次数上限');

end

2.运行结果：

割线方式迭代法

1.程序：

function x=ge_xian_fa(fname,dfname,x0,x1,e,N)

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

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

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

if nargin<5,N=500;end

if nargin<4 e=1e-4;end

while abs(a-b)>e&k

k=k+1;

x=x1-((x1-x0)/(feval(fname,x1)-feval(fname,x0)))*feval(fname,x1);

if feval(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

if k==N,warning('已达迭代次数上限'); end

2.运行结果：