计算方法上机题目.docx
《计算方法上机题目.docx》由会员分享，可在线阅读，更多相关《计算方法上机题目.docx（24页珍藏版）》请在冰豆网上搜索。
计算方法上机题目
1.计算方法A上机作业
上机练习目的
❑复习和巩固数值计算方法的基本数学模型，全面掌握运用计算
机进行数值计算的具体过程及相关问题。
❑利用计算机语言独立编写、调试数值计算方法程序，培养学生
利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务
•利用计算机语言编写并调试一系列数值方法计算通用程序，并
能正确计算给定题目，掌握调试技能。
•掌握文件使用编程技能，如文件的各类操作，数据格式设计、
通用程序运行过程中文件输入输出运行方式设计等。
•写出上机练习报告。
计算方法A上机题目
1.QR分解方法求解线性方程组。
（第二章）
2.共轭梯度法求解线性方程组。
（第三章）
3.三次样条插值（第四章）
4.四阶龙格－库塔法求解常微分方程的初值问题
程序设计要求
1.程序要求是通用的，在程序设计时要充分考虑哪些变量应该可
变的。
2.程序要求调试通过。
上机报告要求
报告内容包括：
●每种方法的算法原理及程序框图。
●程序使用说明。
●算例计算结果。
2.QR分解法求解线性方程组
计算原理
当是任意给定的非零向量，是任意给定的单位向量，则存在初等反射阵,使得，其中为常数，当取单位向量时，由u确定的矩阵H必定满足，所以在计算过程中取u的值为上述值。
设A是一个阶矩阵且它的列向量线性无关，则利用豪斯霍尔德变换可以把A逐步化为上梯形矩阵，设
具体变换过程如下：
设是m维单位坐标向量。
第1步为把矩阵A的第一列化为，取（或取），根据上式可得，取
其中
令，用左乘得，
程序框图
得到矩阵Q和R
i=n
输入系数矩阵A右端列向量b
否
是
x
计算实习
用豪斯霍尔德变换求方程组，其中
Matlab代码
%使用说明：
%需自己输入矩阵A及b的值
%变量Q,R分别为进行QR分解后的结果
clear
clc
formatlong
load（'A矩阵.mat'）
load（'b矩阵.mat'）
%调用函数qrhs进行QR分解
[Q,R]=qrhs（A）;
[~,n]=size（R）;
fprintf（'您输入的矩阵阶数'）;
n
y=Q'*b;
%回代过程
x（n,1）=y（n,1）/R（n,n）;
fori=n-1:
-1:
1
s=y（i,1）;
forj=i+1:
n
s=s-R（i,j）*x（j,1）;
end
x（i,1）=s/R（i,i）;
end
x
被调用函数qrhs
function[Q,R]=qrhs（A）
formatlong
[~,n]=size（A）;
Q=eye（n）;
forj=1:
n-1
B=norm（A（j:
n,j））;
Y=zeros（n-j+1,1）;
Y（1,1）=-sign（A（1,j））*B;
X=A（j:
n,j）;
I=eye（n-j+1）;
N=I-（2/（norm（X-Y））^2）*（X-Y）*（X-Y）';
H=[eye（j-1）zeros（j-1,n-j+1）;zeros（n-j+1,j-1）N];
A=H*A;
Q=Q*H;
end
R=A;
Q;
R;
end
3.
共轭梯度法求解线性方程组
计算原理
当A是n阶对称正定矩阵，若是二次函数的极小值点，则是方程组的解，即
共轭梯度法在形式上具有迭代法的特征，即给定初始向量，由迭代格式
产生的迭代序列在无舍入误差的假定下，最多经过n次迭代就能求得的最小点，也就是方程组的解。
共轭梯度法中关键的两点是，确定迭代格式中的搜索方向和最佳步长。
搜索方向，与前一次的搜索方向关于关于矩阵A共轭，即，然后从点出发，沿方向求得的极小值点，即。
由此解得最佳步长和参数的表达式为
共轭梯度法的计算公式为：
程序框图
计算实习
用共轭梯度法求解线性方程组，其中
矩阵A的阶数n分别取为100,200,400，指出计算结果是否可靠。
Matlab代码
%使用说明：
共轭梯度法求解Ax=b
%命令行中输入矩阵A及b
%然后调用函数getd进行计算
%变量含义：
n—方程阶数，x0—初始向量
%e—计算精度，r—残向量
clear
clc
formatlong
n=input（'请输入方程阶数n='）;%输入矩阵阶数
A=zeros（n）;
b=zeros（n,1）;
A（1,1）=-2;A（1,2）=1;A（n,n-1）=1;A（n,n）=-2;
b（1,1）=-1;b（n,1）=-1;
fori=2:
n-1;
A（i,i）=-2;A（i,i-1）=1;A（i,i+1）=1;
end;
A;
b;%生成对应阶数的矩阵A和b
x0=zeros（n,1）;%生成初始向量
e=input（'请输入计算精度e='）;%输入计算精度
x=getd（A,b,x0,e）;%调用函数getd进行计算
ifn==100%对x元素进行重新排列
x=reshape（x,10,10）
elseifn==200
x=reshape（x,10,20）
else
x=reshape（x,20,20）
end
被调用函数getd
functionx=getd（A,b,x0,e）%矩阵A,b，初始向量x0，精度e
n=size（A,1）;%获取矩阵A的阶数
x=x0;%初始向量
r=b-A*x;%残向量
d=r;%搜索向量
fork=0:
（n-1）
p=（r'*r）/（d'*A*d）;
x=x+p*d;
r2=b-A*x;
if（（norm（r2）<=e）||（k==n-1））
x;
break;
end
q=norm（r2）^2/norm（r）^2;
d=r2+q*d;
r=r2;
end
4.三次样条插值
计算原理
程序框图
计算实习
给定函数.取等距节点，构造三次样条插值。
Matlab代码
%使用说明：
该程序解决的是三次样条插值中，第1,2类边界条件的问题
%各变量含义：
a,b—插值区间左右端点
%n—插值节点数目
%p,q—左右端点导数值
%A，M，d—用于求解AM=d中，矩阵M的值
%C—存放各区间内插值函数的系数矩阵
%zglu—利用追赶法进行LU分解，求解AM=d的函数
clear
clc
formatlong
%输入区间，计算插值节点
a=input（'请输入区间左端点a='）;
b=input（'请输入区间右端点b='）;
n=input（'请输入插值节点数目（包括左右端点）n='）;
d=zeros（n,1）;x=zeros（n,1）;y=zeros（n,1）;A=zeros（n）;
h=（b-a）/（n-1）;
fprintf（'步长h=%d',h）
fori=1:
n
x（i,:
）=a+h*（i-1）;
y（i,:
）=1/（1+25*（x（i,:
））^2）;%计算插值节点处的函数值
end
%选择边界条件进行计算,并输入区间左右端点的导数值p，q
xz=input（'请选择边界条件类型（1或2或3）xz='）;
fprintf（'以第%d类边界条件进行计算',xz）;
p=input（'请输入区间左端边界条件p='）;
q=input（'请输入区间右端边界条件q='）;
%计算矩阵A及矩阵d
ifxz==1
A（1,1）=2;A（1,2）=1/2;A（n,n）=2;A（n,n-1）=1/2;
forj=1:
n-1
d（j,:
）=（3/h）*（（y（j+1,:
）-y（j,:
））/h-（y（j,:
）-y（j-1,:
））/h）;
A（j,j）=2;A（j,j-1）=0.5;A（j,j+1）=0.5;
end
d（1,:
）=d（1,:
）-1/2*p;
d（n,:
）=d（n,:
）-1/2*q
elseifxz==2
d（1,:
）=（6/h）*（（y（2,:
）-y（1,:
））/h-p）;
d（n,:
）=（6/h）*（q-（y（n,:
）-y（n-1,:
））/h）;
A（1,1）=2;A（1,2）=1;A（n,n）=2;A（n,n-1）=1;
forj=2:
n-1
d（j,:
）=（3/h）*（（y（j+1,:
）-y（j,:
））/h-（y（j,:
）-y（j-1,:
））/h）;
A（j,j）=2;A（j,j-1）=0.5;A（j,j+1）=0.5;
end
end
%调用函数zglu用追赶法计算AM=d
M=zglu（A,d）;
%各插值区间内函数表达式，系数矩阵为n*4阶矩阵C
fork=2:
n
C（k-1,1）=（M（k,:
）-M（k-1,:
））/（6*h）;
C（k-1,2）=（x（k,:
）*M（k-1,:
）-x（k-1,:
）*M（k,:
））/（2*h）;
C（k-1,3）=1/（2*h）*（x（k-1,:
）^2*M（k,:
）-x（k,:
）^2*M（k-1,:
））+1/h*（y（k,:
）-1/6*h^2*M（k,:
）-y（k-1,:
）-1/6*h^2*M（k-1,:
））;
C（k-1,4）=1/（6*h）*（x（k,:
）^3*M（k-1,:
）-x（k-1,:
）^3*M（k,:
））+1/h*（x（k,:
）*（y（k-1,:
）-h^2/6*M（k-1,:
））-x（k-1,:
）*（y（k,:
）-h^2/6*M（k,:
）））;
end
%显示输入数据
disp（'您输入的数据如下：
'）
disp（'插值节点x：
'）
x（:
:
）
disp（'插值节点y：
'）
y（:
:
）
disp（'计算得到矩阵M：
'）
M（:
:
）
%输出插值函数S（x）的表达式
disp（'S（x）的表达式为：
'）
forl=1:
n-1
disp（[num2str（C（l,1））,'x^3+',num2str（C（l,2））,'x^2+',num2str（C（l,3））,'x+',num2str（C（l,4））,'',num2str（x（l,:
））,'≤x≤',num2str（x（l+1,:
））]）;
end
被调用函数zglu
%追赶法求解三对角方程组
functionx=zglu（A,b）
[~,n]=size（A）;
L=eye（n）;U=zeros（n）;y=zeros（n,1）;x=zeros（n,1）;
U（1,1）=A（1,1）;y（1,1）=b（1,1）;
fori=2:
n
l=A（i,i-1）/U（i-1,i-1）;
L（i,i-1）=l;
U（i,i）=A（i,i）-l*A（i-1,i）;
y（i,1）=b（i,1）-l*y（i-1,1）;
U（i-1,i）=A（i-1,i）;
end
x（n,1）=y（n,1）/U（n,n）;
fori=n-1:
-1:
1
s=y（i,1）;
forj=i+1:
n
s=s-U（i,j）*x（j,1）;
end
x（i,1）=s/U（i,i）;
end
end
5.四阶龙格－库塔法求解常微分方程的初值问题
计算原理
程序框图
是
输入待求微分方程、求解的自变量范围、初值以及求解范围内的步长等。
确定求解范围内的取点数n
k=n？
否
求解：
求解并输出：
结束算法
计算实习
用标准4级4阶R-K法求解
取步长h=
升级会员 