您好,欢迎来到爱玩科技网。
搜索
您的当前位置:首页数值分析编程及运行结果(高斯顺序消元法)

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

来源:爱玩科技网
高斯消元法

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).

.

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).

.

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')%上面:达到精度要求就结束程序,输出迭代次数和方程组的解 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&kx0=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&kx=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)).

.

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.运行结果:

.

因篇幅问题不能全部显示,请点此查看更多更全内容

Copyright © 2019- aiwanbo.com 版权所有 赣ICP备2024042808号-3

违法及侵权请联系:TEL:199 18 7713 E-MAIL:2724546146@qq.com

本站由北京市万商天勤律师事务所王兴未律师提供法律服务