计算方法上机题目文档格式.docx
- 文档编号:14874380
- 上传时间:2022-10-25
- 格式:DOCX
- 页数:24
- 大小:1.62MB
计算方法上机题目文档格式.docx
《计算方法上机题目文档格式.docx》由会员分享,可在线阅读,更多相关《计算方法上机题目文档格式.docx(24页珍藏版)》请在冰豆网上搜索。
设A是一个阶矩阵且它的列向量线性无关,则利用豪斯霍尔德变换可以把A逐步化为上梯形矩阵,设
具体变换过程如下:
设是m维单位坐标向量。
第1步为把矩阵A的第一列化为,取(或取),根据上式可得,取
其中
令,用左乘得,
程序框图
得到矩阵Q和R
i=n
输入系数矩阵A右端列向量b
否
是
x
计算实习
用豪斯霍尔德变换求方程组,其中
Matlab代码
%使用说明:
%需自己输入矩阵A及b的值
%变量Q,R分别为进行QR分解后的结果
clear
clc
formatlong
load('
A矩阵.mat'
)
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:
s=s-R(i,j)*x(j,1);
end
x(i,1)=s/R(i,i);
end
被调用函数qrhs
function[Q,R]=qrhs(A)
[~,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;
3.
共轭梯度法求解线性方程组
当A是n阶对称正定矩阵,若是二次函数的极小值点,则是方程组的解,即
共轭梯度法在形式上具有迭代法的特征,即给定初始向量,由迭代格式
产生的迭代序列在无舍入误差的假定下,最多经过n次迭代就能求得的最小点,也就是方程组的解。
共轭梯度法中关键的两点是,确定迭代格式中的搜索方向和最佳步长。
搜索方向,与前一次的搜索方向关于关于矩阵A共轭,即,然后从点出发,沿方向求得的极小值点,即。
由此解得最佳步长和参数的表达式为
共轭梯度法的计算公式为:
计算实习
用共轭梯度法求解线性方程组,其中
矩阵A的阶数n分别取为100,200,400,指出计算结果是否可靠。
共轭梯度法求解Ax=b
%命令行中输入矩阵A及b
%然后调用函数getd进行计算
%变量含义:
n—方程阶数,x0—初始向量
%e—计算精度,r—残向量
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)
被调用函数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;
q=norm(r2)^2/norm(r)^2;
d=r2+q*d;
r=r2;
4.三次样条插值
给定函数.取等距节点,构造三次样条插值。
该程序解决的是三次样条插值中,第1,2类边界条件的问题
%各变量含义:
a,b—插值区间左右端点
%n—插值节点数目
%p,q—左右端点导数值
%A,M,d—用于求解AM=d中,矩阵M的值
%C—存放各区间内插值函数的系数矩阵
%zglu—利用追赶法进行LU分解,求解AM=d的函数
%输入区间,计算插值节点
a=input('
请输入区间左端点a='
b=input('
请输入区间右端点b='
请输入插值节点数目(包括左右端点)n='
d=zeros(n,1);
x=zeros(n,1);
y=zeros(n,1);
h=(b-a)/(n-1);
步长h=%d'
h)
fori=1:
x(i,:
)=a+h*(i-1);
y(i,:
)=1/(1+25*(x(i,:
))^2);
%计算插值节点处的函数值
%选择边界条件进行计算,并输入区间左右端点的导数值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:
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;
d(1,:
)=d(1,:
)-1/2*p;
d(n,:
)=d(n,:
)-1/2*q
elseifxz==2
)=(6/h)*((y(2,:
)-y(1,:
))/h-p);
)=(6/h)*(q-(y(n,:
)-y(n-1,:
forj=2:
end
%调用函数zglu用追赶法计算AM=d
M=zglu(A,d);
%各插值区间内函数表达式,系数矩阵为n*4阶矩阵C
fork=2:
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,:
)^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,:
)));
%显示输入数据
disp('
您输入的数据如下:
'
插值节点x:
x(:
:
插值节点y:
y(:
计算得到矩阵M:
M(:
%输出插值函数S(x)的表达式
S(x)的表达式为:
forl=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,:
))]);
被调用函数zglu
%追赶法求解三对角方程组
functionx=zglu(A,b)
L=eye(n);
U=zeros(n);
U(1,1)=A(1,1);
y(1,1)=b(1,1);
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);
x(n,1)=y(n,1)/U(n,n);
s=s-U(i,j)*x(j,1);
x(i,1)=s/U(i,i);
5.四阶龙格-库塔法求解常微分方程的初值问题
输入待求微分方程、求解的自变量范围、初值以及求解范围内的步长等。
确定求解范围内的取点数n
k=n?
求解:
求解并输出:
结束算法
用标准4级4阶R-K法求解
取步长h=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机 题目