数值分析总结
数值分析总结思维导图
Docs
相关代码的使用和注释
列主元Gauss消元法
%%列主元高斯消元法
function x=Gauss_lzy(A,b)%A为方程组系数矩阵，b为方程组的右侧向量，x为方程组的解
[n,m]=size(A);%%得到矩阵A的行和列的宽度
nb=length(b);%%方程组右侧向量的长度
if n~=m%%如果系数矩阵的行数和方程组右侧向量的长度不相等，错误error('%系数矩阵必须是方的');
end
if m~=nb%%方程的变量数和方程右侧向量的长度不相等，错误error('%b的维数与方程的行数不匹配!');
end
for k=1:n-1%%执行n-1次选主元的过程，就可以选完所有，最后剩下的一个直接处理，也表示列%选主元a_max=0;%%先定义一个最大值for i=k:n%%从当前行开始到最后一行选主元if abs(A(i,k))>a_max%%如果遇到比当前最大值大的直接记录作为主元a_max=abs(A(i,k));r=i;%%同时记录下它的行数endendif a_max<1e-15%%如果记录的主元小于1e-5，错误error('%系数矩阵奇异，无法匹配方程组');end%交换两行if r>k%%如果主元所在的行不是当前行，需要交换左侧和右侧for j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%消元过程for i=k+1:n%%从当前行的下一行开始消元m=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);end
end
%回代过程
if abs(A(n,n))<1e-15error('%系数矩阵奇异，无法求解方程组');
end
x=zeros(size(b));
for k=n:-1:1%%从最后一行开始回代for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);
end
Jacobi迭代法
%%Jacobi迭代法
function x=Jacobi(A,b,x0,eps)
D=diag(diag(A));%%将矩阵A的对角元素提取出来
D=inv(D);%%转置
L=tril(A,-1);%%提取矩阵A的下三角
U=triu(A,1);%%提取矩阵A的上三角
B=-D*(L+U);%%雅可比迭代公式
f=D*b;
k=0;
x0=x0;
x=B*x0+f;%%迭代公式
fprintf('k x1_(k) x2_(k) x3_(k)\n');
fprintf('%2d %4.0f %4.0f %4.0f\n',k,x0);
while norm(x-x0)>=eps%%没有到达指定的误差值之前执行循环，不断迭代x0=x;x=B*x0+f;k=k+1;fprintf('%2d %4.0f %4.0f %4.0f\n',k,x0);
endGauss-Seidel迭代法
%%Gauss-seidel迭代法
function x=GaussSeidel(A,b,x0)%%和Jacobi迭代同样的思路，只是公式发生了变化
D=diag(diag(A));
L=tril(A,-1);
C=inv(D+L);
U=triu(A,1);
B=-C*U;
f=C*b;
i=0;
x0=x0;
x=B*x0+f;
fprintf('k x1_(k) x2_(k) x3_(k)\n');
fprintf('%2d %4.0f %4.0f %4.0f\n',k,x0);
for i=1:10x0=x;x=B*x0+f;fprintf('%2d %4.0f %4.0f %4.0f\n',k,x0);
end
二分法
function x = bisectionMethod(A, b, tol)[n, m] = size(A);nb = length(b);if n ~= merror('系数矩阵必须是方的');endif m ~= nberror('b的维数与方程的行数不匹配!');end% 定义二分法的初始下界和上界lower_bound = -1e6;upper_bound = 1e6;% 设置二分法的最大迭代次数max_iterations = 1000;% 循环执行二分法迭代for k = 1:max_iterationslambda = (lower_bound + upper_bound) / 2; % 计算当前迭代的 lambda 值% 解上界对应的方程组并计算残差x_upper = GaussianElimination(A - lambda * eye(n), b);residual_upper = norm(A * x_upper - lambda * x_upper - b);% 解下界对应的方程组并计算残差x_lower = GaussianElimination(A - lower_bound * eye(n), b);residual_lower = norm(A * x_lower - lower_bound * x_lower - b);% 判断是否满足终止条件if abs(residual_upper - residual_lower) < tolbreak;end% 更新下界和上界if residual_upper > residual_lowerupper_bound = lambda;elselower_bound = lambda;endend% 返回最终二分法得到的解x = x_upper;
endfunction x = GaussianElimination(A, b)[n, m] = size(A);nb = length(b);if n ~= merror('系数矩阵必须是方的');endif m ~= nberror('b的维数与方程的行数不匹配!');end% 高斯消元过程for k = 1:n-1% 选主元a_max = abs(A(k, k));r = k;for i = k:nif abs(A(i, k)) > a_maxa_max = abs(A(i, k));r = i;endendif a_max < 1e-15error('系数矩阵奇异，无法匹配方程组');end% 交换两行if r > ktemp = A(k, :);A(k, :) = A(r, :);A(r, :) = temp;temp = b(k);b(k) = b(r);b(r) = temp;end% 消元过程for i = k+1:nm = A(i, k) / A(k, k);for j = k+1:nA(i, j) = A(i, j) - m * A(k, j);endb(i) = b(i) - m * b(k);endend% 回代过程if abs(A(n, n)) < 1e-15error('系数矩阵奇异，无法求解方程组');endx = zeros(size(b));for k = n:-1:1for j = k+1:nb(k) = b(k) - A(k, j) * x(j);endx(k) = b(k) / A(k, k);end
endNewton法
%%Newton法
function x=Newton(fname,dfname,x0,e,N)
%%fname和dfname分别表示f(x)及其导函数的M函数句柄或内嵌函数表达式
if nargin<5,N=500;
end
if nargin<4,e=1e-4;
end
x=x0;
x0=x+2*e;
k=0;
while abs(x0-x)>e&k<N%%大于误差允许值且没有达到迭代次数，继续迭代k=k+1;x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);%%Newton公式fprintf('It.no=%2d x[%2d]=%12.9f\n',k,k,x);
end
if k==N,fprintf('已经达到迭代次数');
endLagrange插值
%%Lagrange插值
function yy=Lagrange(x,y,xi)
m=length(x);%%自变量的长度
n=length(x);%%因变量的长度
if m~=nerror('向量x与y的长度必须一致');
end
s=0;
for i=1:nz=ones(1,length(xi));%%建立一个预备数组for j=1:nif j~=iz=z.*(xi-x(j))/(x(i)-x(j));//%%Lagrange插值公式end
end
s=s+z*y(i);
end
yy=s;