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

高斯消元法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)<e
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)<e
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')<w
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<N,
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<N,
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.运行结果：
Newton插值1.程序：
%保存文件名为New_Int.m
%Newton基本插值公式
%x为向量，全部的插值节点
%y为向量，差值节点处的函数值
%xi为标量，是自变量
%yi为xi出的函数估计值
function yi=newton_chazhi(x,y,xi)
n=length(x);
m=length(y);
if n~=m
error('The lengths of X ang Y must be equal!'); return;
end
%计算均差表Y
Y=zeros(n);
Y(:,1)=y';
for k=1:n-1
for i=1:n-k
if abs(x(i+k)-x(i))<eps
error('the DATA is error!');
return;
end
Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i)); end
end
%计算牛顿插值公式
yi=0;
for i=1:n
z=1;
for k=1:i-1
z=z*(xi-x(k));
end
yi=yi+Y(1,i)*z;
end
2.运行结果：
Lagrange插值1.程序：
function y0 = Language(x,y,x0)
syms t l;
if length(x)==length(y)
n = length(x);
else
disp('x和y的维数不相等！');
return; %检错
end
h=sym(0);
for i=1:n
l=sym(y(i));
for j=1:i-1
l=l*(t-x(j))/(x(i)-x(j));
end;
for j=i+1:n
l=l*(t-x(j))/(x(i)-x(j));
end;
h=h+l;
end
simplify(h);
if nargin == 3
y0 = subs (h,'t',x0); %计算插值点的函数值
else
y0=collect(h);
y0 = vpa(y0,6); %将插值多项式的系数化成6位精度的小数end
2.运行结果：
最小二乘法
1.程序：
function p=nafit(x,y,m)
%多项式拟合
% x, y为已知数据点向量, 分别表示横,纵坐标, m为拟合多项式的次数, 结果返回m次拟合多项式系数, 从高次到低次存放在向量p 中.
A=zeros(m+1,m+1);
for i=0:m
for j=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,'p k');
hold on
plot(t,S3,'r');
2.运行结果：。