高斯消去、追赶法 matlab

1. 分别用Gauss消去法、列主元Gauss消去法、三角分解方法求解方程组

程序:

1Guess消去法:

function x=GaussXQByOrder(A,b)

%Gauss消去法

N = size(A);

n = N(1);

x = zeros(n,1);

for i=1:(n-1)

for j=(i+1):n

if(A(i,i)==0)

disp(‘对角元不能为0‘);

return;

end

m = A(j,i)/A(i,i);

A(j,i:n)=A(j,i:n)-m*A(i,i:n);

b(j)=b(j)-m*b(i);

end

end

x(n)=b(n)/A(n,n);

for i=n-1:-1:1

x(i)=(b(i)-sum(A(i,i+1:n)*x(i+1:n)))/A(i,i);

end

命令行输入:

A=[1 -1 2 1;-1 3 0 -3 ;2 0 9 -6;1 -3 -6 19];

b=[1 3 5 7];

b=b‘;

x=GaussXQByOrder(A,b)

运算结果:

x =

-8.000000000000000

0.333333333333333

3.666666666666667

2.000000000000000

(2)列主元Gauss消去法

程序:

function x=GaussXQLineMain(A,b)

%列主元Gauss消去法

N = size(A);

n = N(1);

x = zeros(n,1);

zz=zeros(1,n);

for i=1:(n-1)

[~,p]=max(abs(A(i:n,i)));

zz=A(i,:);

A(i,:)=A(p+i-1,:);

A(p+i-1,:)=zz;

temp=b(i);

b(i)=b(i+p-1);

b(i+p-1)=temp;

for j=(i+1):n

m = A(j,i)/A(i,i);

A(j,i:n)=A(j,i:n)-m*A(i,i:n);

b(j)=b(j)-m*b(i);

end

end

x(n)=b(n)/A(n,n);

for i=n-1:-1:1

x(i)=(b(i)-sum(A(i,i+1:n)*x(i+1:n)))/A(i,i);

end

命令行:

A=[1 -1 2 1;-1 3 0 -3 ;2 0 9 -6;1 -3 -6 19];

b=[1 3 5 7];

b=b‘;

x=GaussXQLineMain(A,b)

运行结果:

x =

-8.000000000000005

0.333333333333332

3.666666666666668

2.000000000000000

(3)三角分解方法

程序:

function x = LU(A,b)

%三角分解

N = size(A);

n = N(1);

L = eye(n,n);

U = zeros(n,n);

x = zeros(n,1);

y = zeros(n,1);

U(1,1:n) = A(1,1:n);

L(1:n,1) = A(1:n,1)/U(1,1);

for k=2:n

for i=k:n

U(k,i) = A(k,i)-L(k,1:(k-1))*U(1:(k-1),i);

end

for j=(k+1):n

L(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k);

end

end

y(1)=b(1)/L(1,1);

for i=2:n

y(i)=b(i)-sum(L(i,1:i-1)*y(1:i-1));

end

x(n)=y(n)/U(n,n);

for i=n-1:-1:1

x(i)=(y(i)-sum(U(i,i+1:n)*x(i+1:n)))/U(i,i);

end

命令行:

A=[1 -1 2 1;-1 3 0 -3 ;2 0 9 -6;1 -3 -6 19];

b=[1 3 5 7];

b=b‘;

x=LU(A,b)

运行结果:

x =

-8.000000000000000

0.333333333333333

3.666666666666667

2.000000000000000

程序:function [times,wucha]=zhuiganfa(a,b,c,f)

%追赶法:x为所求解,times为所有乘除运算次数(即时间),wucha为误差的2-范数。

n=length(f);

x=zeros(n,1);

y=zeros(n,1);

times=0;

alpha=zeros(1,n);

p=zeros(1,n-1);

alpha(1)=b(1);

for i=2:n

p(i-1)=c(i-1)/alpha(i-1);

alpha(i)=b(i)-a(i-1)*p(i-1);

times=times+1;

end

y(1)=f(1)/b(1);

for i=2:n

y(i)=(f(i)-a(i-1)*y(i-1))/alpha(i);

times=times+1;

end

x(n)=y(n);

for i=n-1:-1:1

x(i)=y(i)-p(i)*x(i+1);

times=times+1;

end

A=zeros(n,n);

A=diag(b,0)+diag(a,-1)+diag(c,1);

wucha=norm((A*x-f‘),2);

命令行(n=20):

a=repmat(11,1,19);

b=repmat(-19,1,20);

c=repmat(7,1,19);

f1=repmat(1.1,1,18);

f=[0 f1 1];

[times,wucha]=zhuiganfa(a,b,c,f)

运行结果:

times =

57

wucha =

8.009010697694412e-15

n=50

命令行:

a=repmat(11,1,49);

b=repmat(-19,1,50);

c=repmat(7,1,49);

f1=repmat(1.1,1,48);

f=[0 f1 1];

[times,wucha]=zhuiganfa(a,b,c,f)

运行结果:

times =

147

wucha =

1.292635294609912e-14

命令行(n=100)

a=repmat(11,1,99);

b=repmat(-19,1,100);

c=repmat(7,1,99);

f1=repmat(1.1,1,98);

f=[0 f1 1];

[times,wucha]=zhuiganfa(a,b,c,f)

结果:

times =

297

wucha =

2.599344850768740e-14

程序:function [count,wucha] = zhouqisanduijaiozhuiganfa(a,b,c,f)

%x为所求解,count为所有乘除运算次数

n=length(f);

x=zeros(n,1);

y=zeros(n,1);

count=0;

l=zeros(1,n-2);

s=zeros(1,n-1);

u=zeros(1,n);

t=zeros(1,n-1);

u(1)=b(1);t(1)=1;

s(1)=1/u(1);y(1)=f(1);

for i=2:n-1

l(i-1)=a(i-1)/u(i-1);

u(i)=b(i)-l(i-1)*c(i-1);

t(i)=-l(i-1)*t(i-1);

s(i)=-s(i-1)*c(i-1)/u(i);

y(i)=f(i)-l(i-1)*y(i-1);

count=count+4;

end

st=0;

for k=1:n-1

st=st+s(k)*t(k);

count=count+1;

end

sy=0;

for k=1:n-2

sy=sy+s(k)*y(k);

count=count+1;

end

u(n)=b(n)-st-s(n-1)*(c(n-1)+t(n-1));

y(n)=f(n)-sy;

x(n)=y(n)/u(n);

for i=n-1:-1:1

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

count=count+1;

end

A=zeros(n,n);

A=diag(b,0)+diag(a,-1)+diag(c,1);

A(n,1)=1;

A(1,n)=1;

wucha=norm((A*x-f‘),2);

命令行:

n=10;

a=repmat(11,1,n-1);b=repmat(-19,1,n);

c=repmat(7,1,n-1);f1=repmat(1.1,1,n-2);f=[0 f1 1];

[count,wucha]= zhouqisanduijaiozhuiganfa(a,b,c,f)

运行结果:

count =

58

wucha =

4.525439045433075

n=30

count =

198

wucha =

5.951269557941316

n=100

count =

688

wucha =

5.993271932634396

原文地址:https://www.cnblogs.com/wander-clouds/p/9991527.html

时间: 2024-11-02 18:18:43

高斯消去、追赶法 matlab的相关文章

计算方法B_高斯消去

%计算方法No.1 %20180916 by wupenghao %高斯消去 %!!!循环中的步长一定要设置准确,+1和-1等,一定要注意!!! A=rand(10,10); b=rand(10,1); x=A\b; %消元 cof=zeros(10,1); root=zeros(10,1); for k=1:1:9 for i=k+1:length(A) cof(i)=A(i,k)/A(k,k); b(i)=b(k)*(-cof(i))+b(i); for j=k:length(A) A(i,

计算方法B_列主元高斯消去

%列主元高斯消去法 %by wu penghao A=rand(10,10); b=rand(10,1); x_c=A\b; %真实值 x=zeros(10,1); n=length(A); %消去过程 for k=1:1:n-1 max=abs(A(k,k)); m=k; for i=k:1:n if max<abs(A(i,k)) max = abs(A(i,k));%每列的最大值 m=i;%每一列最大值索引 end end A([k,m],:)=A([m,k],:);%交换行 b([k,m

图像处理之基础---滤波器 高斯滤波

引用 keendawn 的 高斯(核)函数简介 1函数的基本概念 所谓径向基函数 (Radial Basis Function 简称 RBF), 就是某种沿径向对称的标量函数. 通常定义为空间中任一点x到某一中心xc之间欧氏距离的单调函数 , 可记作 k(||x-xc||), 其作用往往是局部的 , 即当x远离xc时函数取值很小.最常用的径向基函数是高斯核函数 ,形式为 k(||x-xc||)=exp{- ||x-xc||^2/(2*σ)^2) } 其中xc为核函数中心,σ为函数的宽度参数 ,

在实现了高斯消除

高斯消去的实现中使用的增广矩阵成上三角矩阵.然后从下迭代值. 详细是这里. 比方说,有一个线性方程组 然后.出来弄成一个2*2的矩阵,然后再把方程组中等号右边的常数项加进来.成为一个2*3的矩阵 这就是一个增广矩阵了. 接下来变成一个上三角矩阵, 从矩阵的第一行開始,一直到最后一行. 例如说如今面临的是第i行, 那么在i到最后一行,找到第i列数的绝对值最大的那行跟第i行换一个位置,这样交换是有优点的,就是当面对0的情况.哈哈.想一想 然后就用第i行開始把从i+1到最后一行的全部行进行消元操作.消

加性高斯白噪声 AWGN

加性高斯白噪声 AWGN(Additive White Gaussian Noise) 是最基本的噪声与干扰模型. 加性噪声:叠加在信号上的一种噪声,通常记为n(t),而且无论有无信号,噪声n(t)都是始终存在的.因此通常称它为加性噪声或者加性干扰. 白噪声:噪声的功率谱密度在所有的频率上均为一常数,则称这样的噪声为白噪声.如果白噪声取值的概率分布服从高斯分布,则称这样的噪声为高斯白噪声. Matlab中实现加性高斯白噪声: y = awgn(x,SNR) 在信号x中加入高斯白噪声.信噪比SNR

三次样条插值

条件(1)输入$x_{i},y_{i}=f(x_{i}),0\leq i\leq n$(2)要求拟合的曲线$S(x)$满足:对于任意的$1\leq i\leq n-1$,在$x_{i}$处一阶二阶导数连续,$S(x)$ 也连续,且$S^{'}(x_{0})=f^{'}(x_{0})$,$S^{'}(x_{n})=f^{'}(x_{n})$ 求解过程设$S_{x_{j}}^{''}=M_{j}$.对于区间$[x_{j},x_{j+1}]$,$S(x)$是$[x_{j},x_{j+1}]$上的线性函

【原创】开源Math.NET基础数学类库使用(16)C#计算矩阵秩

               本博客所有文章分类的总目录:http://www.cnblogs.com/asxinyu/p/4288836.html 开源Math.NET基础数学类库使用总目录:http://www.cnblogs.com/asxinyu/p/4329737.html 上个月对Math.NET的基本使用进行了介绍,主要内容有矩阵,向量的相关操作,解析数据格式,数值积分,数据统计,相关函数,求解线性方程组以及随机数发生器的相关内容.这个月接着深入发掘Math.NET的各种功能,并对

线性方程组 II

当方程组的未知数个数不等于方程个数时,用高斯消元法得到的是行阶梯型矩阵.此时每个主元所在的列可作为方程组的基本列,基本列的个数为矩阵的秩.选择的列可以不同,但个数唯一.即:当用高斯约当法消减时,可看出非基本列是基本列的线性组合:事实上对线性方程组或者说矩阵的理解有这么几个角度:  1.从行的方向来看 每一行的方程就代表一条直线,解方程组就是找到这些直线的交点. 2.从列的方向来看 可看作列的线性组合,第一列对应第一个未知数,以此类推.方程组的解就是找到这样一组系数,使得矩阵每一列乘以对应未知数系

开源Math.NET基础数学类库使用(16)C#计算矩阵秩

原文:[原创]开源Math.NET基础数学类库使用(16)C#计算矩阵秩                本博客所有文章分类的总目录:http://www.cnblogs.com/asxinyu/p/4288836.html 开源Math.NET基础数学类库使用总目录:http://www.cnblogs.com/asxinyu/p/4329737.html 上个月对Math.NET的基本使用进行了介绍,主要内容有矩阵,向量的相关操作,解析数据格式,数值积分,数据统计,相关函数,求解线性方程组以及