clear all;clc;
fprintf('输入n:(10,20,30)\n');
n=input('');
n
a=zeros(1,n);b=zeros(1,n);c=zeros(1,n-1);
A=zeros(n,n);
for i=1:n
for j=1:n
if j==i
A(i,j)=2;
elseif (j-i==-1)||(j-i==1)
A(i,j)=1;
else
A(i,j)=0;
end
end
end
fprintf('系数矩阵:');A
fprintf('主对角线元素:');b=diag(A)'
fprintf('第1条对角线元素:');c=diag(A,1)'
fprintf('第-1条对角线元素:');a(2:n)=c
d=zeros(1,n);d(1)=1;d(n)=((-1)^(n+1));
fprintf('给定的系数矩阵(转置后):');d
y=zeros(n,1);x=y;
u=zeros(1,n);l=u;%u矩阵U的主对角线元素,l矩阵L的第-1条对角线元素
%---------追过程--------
u(1)=b(1);y(1)=d(1);
for i=2:n
l(i)=a(i)/u(i-1);
u(i)=b(i)-l(i)*c(i-1);
y(i)=d(i)-l(i)*y(i-1);
end
fprintf('中间解向量:');y
fprintf('单位下三角阵L的第-1条对角线元素:');l
fprintf('上三角阵U的主对角线元素:');u
U=zeros(n);
L=eye(n);
for i=1:n-1
L(i+1,i)=l(i+1);
end
fprintf('单位下三角阵:');L %单位下三角阵
for i=1:n-1
U(i,i)=u(i);
U(i,i+1)=c(i);
end
U(n,n)=u(n);
fprintf('上三角阵:');U %上三角阵
%---------赶过程--------
x(n)=y(n)/u(n);
for i=n-1:-1:1
x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
fprintf('解向量:');x
---------------------
作者:zhangchao3322218
来源:CSDN
原文:https://blog.csdn.net/zhangchao3322218/article/details/7412681
版权声明:本文为博主原创文章,转载请附上博文链接!
%追赶法求解三对角线性方程组,Ax=b,A用一维数组a,c,d存储。
function [L,U,x]=crout(a,c,d,b)%数组a存储三角矩阵A的主对角线元素,c、d存储主对角线上边下边带宽为1的元素
n=length(a);
n1=length(c);
n2=length(d);
%错误检查
if n1~=n2%存储矩阵的数组维数错误
error('MATLAB:Crout:不是三对角矩阵,参数数组中元素个数错误.');
elseif n~=n1+1
error('MATLAB:Crout:不是三对角矩阵,参数数组中元素个数错误.');
end
%初始化
L=zeros(n);%生成n*n的全零矩阵
U=zeros(n);
p=1:n;
q=1:n-1;
x=1:n;
y=1:n;
%追赶法程序主体
p(1)=a(1);
for i=1:n-1
q(i)=c(i)/p(i);
p(i+1)=a(i+1)-d(i)*q(i);%d的下标改为1到n-1
end
%正解y
y(1)=b(1)/p(1);%用x存储y
for i=2:n
y(i)=(b(i)-d(i-1)*y(i-1))/p(i);
end
%倒解x
x(n)=y(n);
for i=(n-1):-1:1
x(i)=y(i)-q(i)*x(i+1);
end
%L,U矩阵
for i=1:n
L(i,i)=p(i);
U(i,i)=1;
end
for i=1:n-1
L(i+1,i)=d(i);
U(i,i+1)=q(i);
end %end of function
---------------------
作者:love密密
来源:CSDN
原文:https://blog.csdn.net/jingmiaa/article/details/49660795
版权声明:本文为博主原创文章,转载请附上博文链接!