%
%
%
%% 辣鸡第八题,毁我青春
%% 首先算精确解

clc,clear;
A_matrix=[0.5 1 2;
          1/3 1/4 1/5;
          1/4 1/5 1/6];
B_matrix=[7 8 9]';
RealAns=A_matrix^(-1)*B_matrix
sym(RealAns)

%% 高斯消元法,每一步截取五位有效数字

digits(5) 
Ori_matrix=[0.25 0.2 0.16667 9;
            0.33333 0.25 0.2 8;
            0.5 1 2 7];
vpa(Ori_matrix,21)
FirstStep=double(vpa(OneStepGaoSi(Ori_matrix),5));
vpa(FirstStep,9)%可以看到按照这种操作,可以截取矩阵每个元素保留4位有效数字
FirstStep(2:end,2:end)=double(vpa(OneStepGaoSi(FirstStep(2:end,2:end)),5));
vpa(FirstStep,9)
A=FirstStep(:,1:end-1);
b=FirstStep(:,end);
digits(5) 
x3=vpa(b(3)/A(3,3),5)%用vpa截取x3的5位有效数字
x3=vpa(x3,9)%可以看出直接截取sym型的x3是不行的
x3=vpa(double(x3),5)%先变成double()然后vpa可以截取
vpa(x3,9)%可以看出确实截取了
x2=vpa((b(2)-A(2,3)*x3)/A(2,2),5)
x2=vpa(double(x2),5);
x1=vpa((b(1)-A(1,2)*x2-A(1,3)*x3)/A(1,1),5)
x1=vpa(double(x1),5);
Method1Ans=[x1;x2;x3];

% 列主元消去法

Ori_matrix=[0.5 1 2 7;
            0.33333 0.25 0.2 8;
            0.25 0.2 0.16667 9];
vpa(Ori_matrix,5)
FirstStep=double(vpa(OneStepGaoSi(Ori_matrix),5));
vpa(FirstStep,9)%可以看到按照这种操作,可以截取矩阵每个元素保留4位有效数字
FirstStep(2:end,2:end)=double(vpa(OneStepGaoSi(FirstStep(2:end,2:end)),5));
vpa(FirstStep,9)
A=FirstStep(:,1:end-1);
b=FirstStep(:,end);
x3=vpa(b(3)/A(3,3),5)%用vpa截取x3的5位有效数字
x1=vpa(double(x1),5);
x2=vpa((b(2)-A(2,3)*x3)/A(2,2),5)
x2=vpa(double(x2),5);
x1=vpa((b(1)-A(1,2)*x2-A(1,3)*x3)/A(1,1),5)
x1=vpa(double(x1),5);
Method2Ans=[x1;x2;x3];

% 比较两种解的精确性


norm(RealAns-vpa(Method1Ans,5))
norm(RealAns-vpa(Method2Ans,5))

% 可以看到高斯消元更好
% 以下不用写

y1=A_matrix*Method1Ans
y2=A_matrix*Method2Ans
norm(y1-[9 8 7]')
norm(y2-[9 8 7]')

%

clc,clear;
sym(9/8+3/8*207/434)
sym(3+4527/3472/3)
A=[3 -1 0 0 0;
 -1 3 -1 0 0;
 0 -1 3 -1 0;
 0 0 -1 3 -1;
 0 0 0 -1 3];
b=[1 0 1 0 0]';
A^(-1)*b
[L,U]=lu(A);
U=sym(U')
L=sym(L')
y=U^(-1)*b
L^(-1)*y

%

clc,clear;
syms l21 l31 l32 d1 d2 d3;
L=[1 0 0;
    l21 1 0;
    l31 l32 1];
D=diag([d1,d2,d3]);
L*D*L'
A=[2 -1 1;
-1 -2 3;
1 3 1];
exp=A-L*D*L'==0;
solution=solve(exp,[l21 l31 l32 d1 d2 d3]);
fileds = fieldnames(solution);

for i=1:length(fileds)
    k = fileds(i);
    key = k{1};
    value = solution.(key);
    eval([key,'=value;']);%重要方法,把字符串转化为变量,执行赋值语句
    fprintf('%s=%s\n',key,value)
end
L=[1 0 0;
    l21 1 0;
    l31 l32 1]
D=diag([d1,d2,d3])
b=[6 5 4]';
y=L^(-1)*b
(L*D*L')^(-1)*b

%

clc,clear;
A=[0.6 0.5;0.4 0.3];
ATA=A'*A
simplify(sym((0.36+0.25+0.16+0.09)^0.5))
sym(eig(ATA))
syms l;
eq=det(diag([l,l])-ATA)  
answ=solve(eq,l)
max(answ)^0.5
simplify(max(answ)^0.5)%用simplify分母有理化
vpa(max(answ)^0.5,9)

%

sym([100 99;99 97]^(-1))
A=[100 99;99 97];
ATA=A'*A
syms l;
eq=det(diag([l,l])-ATA);  
answ=solve(eq,l)
double(answ)
(max(answ)/min(answ))^0.5
simplify((max(answ)/min(answ))^0.5)
vpa(double((max(answ)/min(answ))^0.5),9)
syms l;
eq=det(diag([l,l])-A);  
answ=solve(eq,l)
double(answ)
(max(answ)/min(answ))
simplify((max(answ)/min(answ)))

%

simplify(max(answ)/min(answ))
syms a b;
ep = (a*sqrt(4357)+b)^2==simplify(max(answ)/min(answ))
ab=solve(ep,[a b])
ab.b
a=simplify((max(answ)/min(answ))^0.5)
a=simplify(a^0.5)
a=simplify(a^0.5)
a=simplify(a^0.5)
expand(a^2)