计算方法上机题目

目录1.计算方法A 上机作业 (1)上机练习目的 (1)上机练习任务 (1)计算方法A 上机题目 (1)程序设计要求 (1)上机报告要求 (1)2.QR 分解法求解线性方程组 (2)计算原理 (2)程序框图 (7)计算实习 (8)Matlab代码 (8)3.共轭梯度法求解线性方程组 (10)计算原理 (10)程序框图 (11)计算实习 (12)Matlab代码 (12)4.三次样条插值 (14)计算原理 (14)程序框图 (16)计算实习 (17)Matlab代码 (17)5.四阶龙格－库塔法求解常微分方程的初值问题 (21)计算原理 (21)程序框图 (22)计算实习 (23)Matlab代码 (23)1.计算方法A 上机作业上机练习目的❑ 复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

❑ 利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。

上机练习任务•利用计算机语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。

•掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。

•写出上机练习报告。

计算方法A 上机题目1. QR 分解方法求解线性方程组。

（第二章）2. 共轭梯度法求解线性方程组。

（第三章）3. 三次样条插值（第四章）4. 四阶龙格－库塔法求解常微分方程的初值问题程序设计要求1. 程序要求是通用的，在程序设计时要充分考虑哪些变量应该可变的。

2. 程序要求调试通过。

上机报告要求报告内容包括：● 每种方法的算法原理及程序框图。

● 程序使用说明。

● 算例计算结果。

2. QR 分解法求解线性方程组计算原理当nx R∈是任意给定的非零向量，nv R∈是任意给定的单位向量，则存在初等反射阵2THI u u=-,使得H xvσ=，其中σ为常数，当取单位向量2x v ux vσσ-=-时，由u 确定的矩阵H 必定满足H xvσ=，所以在计算过程中取u 的值为上述值。

设A 是一个()m n mn ⨯≥阶矩阵且它的列向量线性无关，则利用豪斯霍尔德变换可以把A 逐步化为上梯形矩阵，设()11121212221212,,,n nnm m m n a a a a a a A a a a a a a ⎛⎫⎪⎪== ⎪ ⎪⎝⎭具体变换过程如下：设()1,2,,i e in =是m 维单位坐标向量。

第1步 为把矩阵A 的第一列()01a 化为()1,0,,0Tσ，取()()011,1,0,,0Tmx a v e R===∈ （或取1ve =-），根据上式可得，取()()0111110121112a e u a e σωωσ-==-其中()()001211121111, a e a e ωσωσ=-=-()()()()01111111111120121111122, T TTTu ua a ωωωωωωασσωασσ====--令111111122TTm m H I u u I αωω-=-=-，用1H 左乘()0A 得，()()()()()()()()()()100001111211111123,,, =,,,,n nAH AH a H a H a e a a a σ==程序框图22i i i i i iTa e I u uσ--计算实习用豪斯霍尔德变换求方程组xA b = ，其中54756753941287886537810987756, 579119755368891089587877810105756759101052A b ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦Matlab 代码 %使用说明：%需自己输入矩阵A 及b 的值%变量Q,R 分别为进行QR 分解后的结果 clear clc format long 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); for i=n-1:-1:1 s=y(i,1);for j=i+1:ns=s-R(i,j)*x(j,1);endx(i,1)=s/R(i,i);endx被调用函数qrhsfunction[Q,R]=qrhs(A)format long[~,n]=size(A);Q=eye(n);for j=1:n-1B=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;endR=A;Q;R;end3. 共轭梯度法求解线性方程组计算原理当A 是n 阶对称正定矩阵，若*x 是二次函数()12T Tf x x A x b x=-的极小值点，则*x 是方程组A xb =的解，即()()n**x Rm in A x b fx fx ∈=⇔=共轭梯度法在形式上具有迭代法的特征，即给定初始向量(0)x ，由迭代格式()()()1k k k k xxdα+=+产生的迭代序列在无舍入误差的假定下，最多经过n 次迭代就能求得()f x 的最小点，也就是方程组A xb =的解。

共轭梯度法中关键的两点是，确定迭代格式中的搜索方向和最佳步长。

搜索方向()kd ，与前一次的搜索方向关于()1k d -关于矩阵A 共轭，即()()()1,0k k d A d+=，然后从点()k x出发，沿()k d方向求得()fx 的极小值点()1k x+，即()()()()()()0m in k k kk fx df x dααα≥+=+ 。

由此解得最佳步长k α和参数k β的表达式为()()()()()()()()1, k Tk k T k kkk Tk k Tk rdrA d d A ddA dαβ+==-共轭梯度法的计算公式为：()()()()()()()()()()()()()()()()()00011(1)11(1)r rr k k T k k k T k k k k k k k T k k k T k k k k k d b A x r d d A dx x d b A x rA d d A d dd ααββ++++++⎧==-⎪⎪=⎪⎪⎪=+⎪⎨=-⎪⎪⎪=-⎪⎪=+⎪⎩程序框图计算实习用共轭梯度法求解线性方程组xA b=，其中211121, 1210121A b --⎡⎤⎡⎤⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦矩阵A 的阶数n 分别取为100,200,400，指出计算结果是否可靠。

Matlab 代码%使用说明：共轭梯度法求解Ax=b %命令行中输入矩阵A 及b %然后调用函数getd 进行计算%变量含义：n —方程阶数，x0—初始向量 %e —计算精度，r —残向量 clear clc format longn=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; for i=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 进行计算if n==100 %对x元素进行重新排列x=reshape(x,10,10)elseif n==200x=reshape(x,10,20)elsex=reshape(x,20,20)end被调用函数getdfunction x=getd(A,b,x0,e)% 矩阵A,b，初始向量x0，精度e n=size(A,1);% 获取矩阵A的阶数x=x0;%初始向量r=b-A*x;%残向量d=r;%搜索向量for k=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;endq=norm(r2)^2/norm(r)^2;d=r2+q*d;r=r2;end4.三次样条插值计算原理程序框图计算实习 给定函数21(x )(1x 1)115f x=-≤≤+ .取等距节点，构造三次样条插值10(x )S 。

Matlab 代码%使用说明：该程序解决的是三次样条插值中，第1,2类边界条件的问题 %各变量含义：a,b —插值区间左右端点 % n —插值节点数目 % p,q —左右端点导数值% A ，M ，d —用于求解AM=d 中，矩阵M 的值 % C —存放各区间内插值函数的系数矩阵% zglu —利用追赶法进行LU 分解，求解AM=d 的函数 clear clc format long%输入区间，计算插值节点 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) for i=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及矩阵dif xz==1A(1,1)=2;A(1,2)=1/2;A(n,n)=2;A(n,n-1)=1/2;for j=1:n-1d(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;endd(1,:)=d(1,:)-1/2*p;d(n,:)=d(n,:)-1/2*qelseif xz==2d(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;for j=2:n-1d(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;endend%调用函数zglu用追赶法计算AM=dM=zglu(A,d);%各插值区间内函数表达式，系数矩阵为n*4阶矩阵Cfor k=2:nC(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)的表达式为：')for l=1:n-1disp([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% 追赶法求解三对角方程组function x=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);for i=2:nl=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);endx(n,1)=y(n,1)/U(n,n);for i=n-1:-1:1s=y(i,1);for j=i+1:ns=s-U(i,j)*x(j,1);endx(i,1)=s/U(i,i);endend5.四阶龙格－库塔法求解常微分方程的初值问题计算原理程序框图计算实习用标准4级4阶R-K 法求解23,(0)1,y (0)3,y (0)2y y y y x y ''''''=+-+-⎧⎨'''=-==⎩取步长h=0.05，计算(1)y 的近似值，并与解析解(x )x e 21xy x =+-作比较。