matlab方程组迭代法,matlab 求解线性方程组——迭代法

标签(空格分隔)： matlab 线性方程组 迭代 数值

文章目录标签(空格分隔)： matlab 线性方程组 迭代 数值@[toc]matlab 求解线性方程组——迭代法1 测试迭代法求解函数GaussSeidelIteration()、JacobiIteration()2 迭代法求解函数GaussSeidelIteration()3 迭代法求解函数JacobiIteration()4 求解矩阵的范数matrixnorm()5 求解矩阵的范数条件数matrixcond()6 向量的p范数vectornorm()联系作者 definedone@163.com

matlab 求解线性方程组——迭代法

1 测试迭代法求解函数GaussSeidelIteration()、JacobiIteration()

clc,clear

% m=[1 2

% -1 2]

% n1=matrixnorm([1 2;-1 2],1)

% n2=matrixnorm([1 2;-1 2],2)

% ninf=matrixnorm([1 2;-1 2],inf)

%

% cond1=matrixcond(m,1)

% cond2=matrixcond(m,2)

% condinf=matrixcond(m,inf)

% a=[10.2 -0.25 -.3

% -.41 1.13 -.15

% -.25 -.14 1.21]

% b=[.515

% 1.555

% 2.78]

% precision=1e-3

% x=[2 2 2]'

% x=JacobiIteration(a,b,precision,x)

% xsystem=inv(a)*b

a=[6 -1 -1

-1 6 -1

-1 -1 6]

b=[11.33 32 42]'

precision=1e-3

x=[4 6 8]'

x=GaussSeidelIteration(a,b,precision,x)

xsystem=inv(a)*b

运行结果：

iterationTimes =

7

x =

4.6661

7.6189

9.0475

xsystem =

4.6661

7.6189

9.0475

2 迭代法求解函数GaussSeidelIteration()

function x=GaussSeidelIteration(a,b,precision,x)

%时间：2013/12/08 姓名：邓能财

disp('%GaussSeidel——迭代法求解线性方程组')

if nargin<4, x=zeros(dim,1);

elseif nargin<3, precision=1e-9; end

[dim,dim2]=size(a);

disp('迭代的方程为：x(k+1)=s*xk+f')

%错误的矩阵：

assert( dim==dim2 && dim==length(b),...

['Argument input error: ',...

'Matrix dimensions must agree.'])

d=diag(diag(a))

a=a-d;

l=-tril(a)

u=-triu(a)

invd_l=inv(d-l)

s=invd_l*u

f=invd_l*b

x_=ones(dim,1);

precision=precision/10;

x

iterationTimes=0

while any(abs(x-x_)>=precision)

x_=x;

x=s*x_+f

iterationTimes=iterationTimes+1

end

end

3 迭代法求解函数JacobiIteration()

function x=JacobiIteration(a,b,precision,x)

%Jacobi迭代法求解线性方程组

%时间：2013/12/08 姓名：邓能财

[dim,dim2]=size(a);

if nargin<4, x=zeros(dim,1);

elseif nargin<3, precision=1e-10; end

disp('迭代的方程为：x(k+1)=g*xk+f')

%错误的矩阵：

assert( dim==dim2 && dim==length(b),...

['Argument input error: ',...

'Matrix dimensions must agree.'])

d=diag(diag(a))

invd=inv(d)

g=invd*(d-a)

f=invd*b

x_=ones(dim,1);

iterationTimes=0

while any(abs(x-x_)>=precision)

x_=x;

x=g*x_+f

iterationTimes=iterationTimes+1

end

end

4 求解矩阵的范数matrixnorm()

function n=matrixnorm(m,p)

%求矩阵的p(1,2,inf)范数

%时间：2013/12/08 姓名：邓能财

if nargin<2, p=1; end

[dim,dim2]=size(m);

%错误的矩阵：

assert( dim==dim2,...

['Argument input error: ',...

'Matrix dimensions must agree.'])

switch p

case 1

disp('矩阵的一范数')

sum_=sum(m,1)

ssum_=sym(sum_)

[max_,maxindex]=max(sum_)

n=max_

case inf

disp('矩阵的无穷范数')

sum_=sum(m,2)

ssum_=sym(sum_)

[max_,maxindex]=max(sum_)

n=max_

case 2

disp('矩阵的2范数')

t=m'*m

st=sym(m'*m)

%t的特征多项式

syms x

tx=x*eye(dim)-t

fx=det(tx)

eigenvalue=eig(t)

seig=eig(st)

sym(eigenvalue)

real_eigenvalue=real(abs(eigenvalue))

[max_,maxindex]=max(real_eigenvalue)

n=sqrt(max_)

otherwise

error('第三个参数只能为1、2或inf')

end

end

5 求解矩阵的范数条件数matrixcond()

function c=matrixcond(m,p)

%时间：2013/12/08 姓名：邓能财

if nargin<2, p=1; end

[dim,dim2]=size(m);

%错误的矩阵：

assert( dim==dim2,...

['Argument input error: ',...

'Matrix dimensions must agree.'])

switch p

case 1

disp('矩阵的一范数条件数')

inv_m=inv(m)

sinv_m=inv(sym(m))

c=matrixnorm(m,p)*matrixnorm(inv_m,p)

case inf

disp('矩阵的无穷范数条件数')

inv_m=inv(m)

sinv_m=inv(sym(m))

c=matrixnorm(m,p)*matrixnorm(inv_m,p)

case 2

disp('矩阵的2范数条件数')

inv_m=inv(m)

sinv_m=inv(sym(m))

c=matrixnorm(m,p)*matrixnorm(inv_m,p)

otherwise

error('第二个参数只能为1、2或inf')

end

end

6 向量的p范数vectornorm()

function n=vectornorm(x,p) %向量的p范数

%时间：2013/12/08 姓名：邓能财

if nargin<2, p=1; end

%错误的p：

assert( p>0,...

['Argument input error: ',...

'The next argument must be positive number! '])

n=sum(abs(x).^p)^(1/p);

end

联系作者 definedone@163.com

end