计算方法A上机作业

研究生课程计算方法A上机作业： 1.QR分解 2.共轭梯度法 3.三次样条 4.龙格-库塔

%%%% --------共轭梯度法求解线性方程组---------%%%%
% %% 参数说明
% A   --系数矩阵，方程Ax=b中的A矩阵，人为给定
% b   --方程Ax=b中的b矩阵，人为给定
% x0  --迭代初始值，人为给定
% e   --给定精度ε，人为给定
% n   --迭代次数，人为给定
% r   --残向量r（k）
% r_  --残向量r（k+1）
% a   --搜索步长α
% bt  --代表β
% norm_r   --r（k）的2范数
% norm_r_  --r（k+1）的2范数
close all; clear all; clc
disp（‘-----共轭梯度法求解线性方程组-----‘）
prompt1=‘请输入系数矩阵A\n‘;
A = input（prompt1）;                     % 输入一个矩阵用于计算
%% 本算法中检测了矩阵A是否为对称正定矩阵，矩阵A对称非正定矩阵，-A为对称正定矩阵，本例题计算得为 -Ax = -b;
e = 0.001;
n = 500;
[m_An_A] = size（A）;                   % 获取矩阵的维数
flag = 1;                            % 标志变量，用于检测系数矩阵A是否为对称正定矩阵
%% 判断矩阵A是否是对称正定矩阵
if m_A == n_A
    for q = 1:length（A）
        c = mat2cell（A[qlength（A）-q][qlength（A）-q]）;
        if det（c{11}）<=0
            flag = 0  ;                 % 如果顺序主子式的行列式小于0，flag置0
        end
    end        
else
    flag = 0;
end
%% 计算
if flag == 0                           
    disp（‘错误！矩阵A不是对称正定矩阵。‘）    % 如果A不是对称正定矩阵，不执行迭代操作，输出提示文字
else
    prompt2=‘请输入对应的矩阵b\n‘;
    b = input（prompt2）;                     % 输入一个矩阵b用于计算
    prompt3=‘请给定初始向量x0\n‘;
    x0 = input（prompt3）;                    % 输入一个初始向量x0
    prompt4=‘请给定精度ε\n‘;
    e = input（prompt4）;                     % 输入一个给定精度ε，精度满足要求是跳出循环
    prompt5=‘请给定迭代次数\n‘;
    n = input（prompt5） ;                    % 输入一个对大迭代次数n  
    r_ = b-A*x0;                           % 计算r0，此处r_代表r0
    d = r_;                                % 计算d0，此处d代表d0
    x = x0;                   
 %% 迭代n次
    for k = 0:n-1
        norm_r_=0;                         % 初始化变量，用于计算r（k）和r（k+1）的2范数，注意这个初始化一定要放在迭代循环里边。
        norm_r =0;
        r = r_;
        a = r‘*r /（d‘*A*d）;
        x = x + a*d;
        r_ = b-A*x;
        %% 计算r（k）和r（k+1）的2范数
        for i_r = 1:m_A
            norm_r_ = norm_r_+ r_（i_r1）*r_（i_r1）;
            norm_r  = norm_r + r（i_r1）*r（i_r1）;
        end
        norm_r_ = sqrt（norm_r_）;
        norm_r  = sqrt（norm_r）;
        bt = norm_r_*norm_r_/（norm_r*norm_r）;    
        d = r_ + bt*d;
        %% 判断精度是否达到要求或迭代次数是否到达上限
        if norm_r_<= e || k+1 == n
            break
        end 
    end   
    disp（[‘残差的2范数为：‘num2str（norm_r_）]）
    disp（[‘迭代‘num2str（k+1）‘次求得方程的解为：‘]）
    x
end


%%% P113 3.2计算实习参数
% n=100 A = -2*eye（100）+[zeros（1100）;eye（99）zeros（991）]+[zeros（991）eye（99）;zeros（1100）]
% n=200 A = -2*eye（200）+[zeros（1200）;eye（199）zeros（1991）]+[zeros（1991）eye（199）;zeros（1200）]
% n=400 A = -2*eye（400）+[zeros（1400）;eye（399）zeros（3991）]+[zeros（3991）eye（399）;zeros（1400）]
% n=100 b = [-1;zeros（981）;-1]    
% n=200 b = [-1;zeros（1981）;-1] 
% n=400 b = [-1;zeros（3981）;-1]     

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3124  2019-01-05 23:32  程序\Conjugate_Gradient.m

     文件        700  2019-01-03 14:47  程序\NewtonInterpolation.m

     文件        670  2019-01-03 14:51  程序\NewtonInterpolation1.m

     文件       2245  2019-01-05 23:09  程序\QR.m

     文件       1262  2019-01-06 09:09  程序\RungeKutta.m

     文件        149  2019-01-06 10:19  程序\sf.m

     文件       6032  2019-01-06 08:36  程序\Spline.m

     文件        456  2018-12-26 22:35  程序\ZG.m

     目录          0  2019-01-05 23:35  程序

----------- ---------  ---------- -----  ----

                14638                    9