计算方法上机题目

《计算方法上机题目》由会员分享，可在线阅读，更多相关《计算方法上机题目（32页珍藏版）》请在装配图网上搜索。

1、目录1. 计算方法 A上机作业1.上机练习目的1.上机练习任务1.计算方法 A 上机题目 1.程序设计要求1.上机报告要求1.2. QR分解法求解线性方程组 2计算原理2.程序框图8.计算实习9.Matlab 代码 9.3. 共轭梯度法求解线性方程组11计算原理1.1程序框图12计算实习13Matlab 代码 134. 三次样条插值1.5计算原理15程序框图17计算实习18Matlab 代码 185. 四阶龙格库塔法求解常微分方程的初值问题 2 2计算原理22程序框图24计算实习25Matlab 代码251. 计算方法A上机作业上机练习目的复习和巩固数值计算方法的基本数学模型，全面掌握运用计算

2、 机进行数值计算的具体过程及相关问题。利用计算机语言独立编写、调试数值计算方法程序，培养学生 利用计算机和所学理论知识分析解决实际问题的能力。上机练习任务?利用计算机语言编写并调试一系列数值方法计算通用程序，并 能正确计算给定题目，掌握调试技能。?掌握文件使用编程技能，如文件的各类操作，数据格式设计、 通用程序运行过程中文件输入输出运行方式设计等。?写出上机练习报告。计算方法A上机题目1. QR分解方法求解线性方程组。（第二章）2. 共轭梯度法求解线性方程组。（第三章）3. 三次样条插值（第四章）4. 四阶龙格-库塔法求解常微分方程的初值问题程序设计要求1. 程序要求是通用的，在程序设计时要充

3、分考虑哪些变量应该可 变的。2. 程序要求调试通过。上机报告要求报告内容包括：每种方法的算法原理及程序框图。程序使用说明。算例计算结果。2. QR分解法求解线性方程组计算原理当x Rn是任意给定的非零向量，v Rn是任意给定的单位向量，则存在初等反射阵H I 2uuT，使得Hxv,其中为常数，当取单位向量ux vPx v F2时，由u确定的矩阵H必定满足Hx v，所以在计算过程中取u的值为上述值 设A是一个m n m n阶矩阵且它的列向量线性无关，则利用豪斯霍尔德a11a12La1na?1A21Ma22La2nMMam1am2Lamn具体变换过程如下：变换可以把A逐步化为上梯形矩阵，设设e i

4、 1,2,L ,n是m维单位坐标向量ai , a2 丄,a第1步为把矩阵A的第一列a 化为仆0丄，0 丁 ，取0x ai ,v e1,0,L ,0 T Rm （或取vei），根据上式可得，取其中2uuTU10a11e1Pa101e PP 1Pp 1 P Pa01eP, 10a1 1eT1 1T1 12 T1 1PP21 2011 a101 1 a11歸白誉-ffi 支 SIF 弓厉ofM f 匚 - -B 匚6).屈 J (E; %)佥 % O 1 (VU ? 丿(首摄口 JT USTSuro cHL L Lue T EeH ).构造用Hk左乘4（i】）得占=濟川2丄）=WkHfc_iHXA（

5、*匸 rarenai + 血呱，力 a黒与+ 叶帕此名,GiV础亠i卫alk起+35必-1Kfl(3)a2k叫Hl*(2)a2.fi+2，-*07JJi-#Vi * bSU1J*)afc.k+2awak+.Jt+l(k)aU12 ft申awam.klr其中（4? 4?皿疋护时=（1+4 - 円刊琬T=耳$7 -夠他二站D -為（聲7 叶占7）, j = 2 I,严由此即得i-*+l)继续上述变换过程，直至s = minm - 1加步，最后得到A = HsH9-H1A=R.曲于每个都是正交矩陆所以H =匕比血也是正交矩阵.上式可以写 为HA = R.当m=n时，R是n阶上三角矩阵.这时，式(2,

6、4.6)可以写为A = 其中Q = H1为n阶正交矩阵，即把n阶矩阵A分解为一个n阶正交矩阵与 个可逆的n阶上三角矩阵R的乘积程序框图:输入系数矩阵 A右端列向量by QTbRx y用豪斯霍尔德变换求方程组Axb其中54756753941287886537810987756A 57911975 ,b 53688910895878778101057计算实习5675910 1052Matlab代码%使用说明：%需自己输入矩阵A及b的值%变量Q,R分别为进行QR分解后的结果clearclcformat Io ngload(A 矩阵.mat)load(b 矩阵.mat)%调用函数qrhs进行QR分解Q

7、,R=qrhs(A);卜,n =size(R);fprintf(您输入的矩阵阶数);ny=Q*b;%回代过程x( n,1)=y( n,1)/R( n,n);s=y(i,1);for j=i+1: ns=s-R(i,j)*x(j,1);endx(i,1)=s/R(i,i);endx被调用函数qrhsfun ctio nQ,R=qrhs(A)format Io ng, n=size(A);Q=eye( n);for j=1: n-1B=n orm(A(j: n,j);Y=zeros( n-j+1,1);Y(1,1)=-sig n(A(1,j)*B;X=A(j: n,j);I=eye( n-j+1)

8、;N=l-(2/( norm(X-Y)F2)*(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*是二次函数f x -xTAx bTx的极小值点, 2则x*是方程组Ax b的解，即Ax b f x min f xn共轭梯度法在形式上具有迭代法的特征，即给定初始向量 x(0)，由迭代格式kd产生的迭代序列在无舍入误差的假定下，最多经过n次迭代就能求得f x的最小点，也就是方程组Ax b的解。共轭梯度法中关键的两点

9、是，确定迭代格式中的搜索方向和最佳步长。 搜索方向d k，与前一次的搜索方向关于d k 1关于矩阵A共轭，即d k , Ad k 10，然后从点x k出发，沿d k方向求得fx的极小值点x k 1kkx d min f0由此解得最佳步长k和参数k的表达式为k T , kr dk TkAdk 1 T k r Ad,k T , k d Ad共轭梯度法的计算公式为:d0Ax0kkx(k 1) rdk1br kTr d k T 典.k d Adkkx kdk 1b Axk 1 T , kr Ad,k T kd Ad(k 1)krkd程序框图计算实习用共轭梯度法求解线性方程组Ax b，其中21112 1

10、0A0 0 0 ,bM1 2 101 21矩阵A的阶数n分别取为100,200,400,指出计算结果是否可靠Matlab代码%使用说明：共轭梯度法求解Ax=b%命令行中输入矩阵A及b%然后调用函数getd进行计算%变量含义：n方程阶数，x0初始向量%e 计算精度，r残向量clearclcformat Io ngn=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,

11、i-1)=1;A(i,i+1)=1;end;A;b;%生成对应阶数的矩阵 A和bx0=zeros( n,1);%生成初始向量e=input(请输入计算精度e=);%输入计算精度%调用函数getd进行计算%对x元素进行重新排列x=getd(A,b,xO,e);if n=100 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=

12、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 (n orm( r2) =e)|(k=n-1)x;break;endq=n orm(r2F2/norm(rF2;d=r2+q*d;r=r2;end4. 三次样条插值计算原理设在区间臥b上给定个If点曲oisiii *，在节点斗Jt的函数值为0丄从若函数S(X)満足以下三个条杵I(1)在每个子区何国亠知丄川)上乂”是三次荽项式I 5(X/) = ViO-0在区间0上，SCO的二阶导数亍连续.则称SCr)为函数眇二An在区ffla.d上的三次样籍插tiSft.

13、由定义可feJ5(x)共有钿个持定参数，根据条杵G)可碍如下dm3个方樫.工(马)吐0再60上条件(2)的”1卜方程”共有4n-2个方程，因此还需嬰境加两个条杵才可以确定 出3亍待定罄敷.所増加的条fl称为边界条佯或爛点杀件常用的三种边界条件为* er(fl)y(*)：已rw-r(6)：已知yw是以r-i为周期的周期函设S(Q在节点曲社的二阶导数值为M(i = 0丄刃)何为持毘参数根据公式(5)可 以得到scoa区间砂】对上的表达式为聚得到的的表达氏*就要确定M的惟对于平同的边界条件，求胖M的方程iflrt(7H9)(2 1心4/h 2易AMV胚d2Pt 21心4g *心佥见d!i 1d 如【

14、丿V M1/Mh:,心X陆(9T-* A -1-M*/-ll -.H-l.公式(8)中的彳丸/Ij.xJ公式中的如一心=心片7儿-Ft三伙样条插值嗜数的关犍是要解式(7). (8)JS(9所示的方程组.椿帯到的JWftA 公武(中可稱和相战的插值霸截.程序框图计算实习1给定函数f(x)- 2( 1 x 1) 取等距节点，构造三次样条插值 S!0(x)115xMatlab代码%使用说明：该程序解决的是三次样条插值中，第1,2类边界条件的问题%各变量含义：a,b插值区间左右端点%n插值节点数目%p,q左右端点导数值%A，M，d用于求解AM=d中，矩阵M的值%C存放各区间内插值函数的系数矩阵%zgl

15、u利用追赶法进行LU分解，求解AM=d的函数 clearclcformat lo ng%输入区间，计算插值节点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: nx(i,:)= a+h*(i-1);y(i,:)=1/(1+25*(x(i,:)A2);%计算插值节点处的函数值end%选择边界条件进行计算，并输入区间左右

16、端点的导数值p，qxz=input(请选择边界条件类型(1或2或3)xz=);fprintf(以第d类边界条件进行计算：xz);p=i nput(请输入区间左端边界条件p=);q=i nput(请输入区间右端边界条件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

17、,:)=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阶

18、矩阵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,:)A2*M(k,:)-x (k, :)A2*M(k-1,:)+1/h*(y(k,:)-1/6*hA2* M(k,:)- y(k-1,:)-1/6*hA2*M(k-1,:);C(k-1,4)=1/(6*h)*(x (k,:)A3*M(k-1,:)-x(k-1,:)A3*M(k,:)+1/h*(x (k, :)*(y(k-1,:)-hA2/6*M(k-1,:)-

19、x(k-1,:)*(y (k,:)-hA2/6* M(k,:);end%显示输入数据disp(您输入的数据如下：)disp(插值节点x：)x(:,:)disp(插值节点y：)y(:,:)disp(计算得到矩阵M :)M(:,:)%输出插值函数S(x)的表达式disp(S(x)的表达式为：)for 1=1: n-1disp( n um2str(C(l,1),xA3+ ,n um2str(C(l,2),xA2+, nu m2str(C(l,3),x+ ,n um2str(C(l,4),num2str(x(l,:), x ,num2str(x(l+1,:);end被调用函数zglu%追赶法求解三对角

20、方程组fun cti on 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

21、(j,1);endx(i,1)=s/U(i,i);endend5. 四阶龙格-库塔法求解常微分方程的初值问题计算原理利用泰勒展开可以导出龙格-库塔（fttinge-Kutta）法（I称R-K方法），m级 的龙格-库塔法的一般形式为Vil =恥 + 入+ 入2K7 + + 人mKuMK =4 Kg - hf（ii 也处 +/?21 Kh）tKr = hf （J； + &i Ki + /?322）,L Km = hfXi + Ctrnhi 班 + ZTTll A 1 + EtTIKh T - + ZSmbTTl-lF71- 1),其中UzB讣均为常数.由待定系数法确定.确定的原则是将局瑚截斷谋差/

22、?# = 讥欧+J-册在航处泰勒展开，适当选取h的系数，使得局部截断误差Ry的阶 尽可能高一类似于匕述二阶方法的推导，可得多种4级4阶FbK法*应用最广泛的是如 下标准（经典）4级4阶RK海彷七=站+孑(皿+ 2也十2/Cj + Ki)t DKi = hfztyi),(9.L20)K? = hf (野+瓠葩+詛)f= hf(xi +h,vt+ 岭),E阶常微分方程初值问趣为円)=f(工必环】小盘冬冬久(0 =州如=如，严7时=牝小艮求解高阶微分方程初值问题是将其转化为阶微分方程组来求解.为此-引进新的 变量班=如宓=讥、加=射-即可将m阶微分方程转化为如下 的一阶微分方程组：V =临(U.3.

23、6)城=阪i/Jn-1 = VmjVm = f(斗如厂伽h.Vi(a)=加血工気加(a) 疔1 1J-程序框图计算实习用标准4级4阶R-K法求解y y y y 2x 3,y(0)1,y(0)3,y(0)2取步长h=0.05,计算y(1)的近似值，并与解析解y(x) xex 2x 1作比较Matlab代码%使用说明：%变量含义：m 微分方程阶数 a, b计算的区间端点%h步长%fm 微分方程的表达形式，按以下形式输入y(l,l)=y,y(2,1)=y,y(3,l)=y clearclc format Io ngm=i nput(请输入常微分方程的阶数m=);a=input(请输入x下限a=);b

24、=input(请输入x上限b=);h=input(请输入步长h=);fm=input(令 y(1,1)=y,y(2,1)=y ,y(3,1)=y 请输入 ym=,s);%微分方程按以下形式输入y(1,1)=y,y(2,1)=y,y(3,1)=y%以下计算过程为一阶初值问题if m=1mm=(b-a)/h;y(1,1)=i nput(请输入在初值点的函数值f(a)=);x=a;y11(1)=y(1,1);for k1=2:(mm+1)%计算K1y1=y(1,1);K(1,1)=h*(eval(fm); x=x+h/2;y(1,1)=y1+K(1,1)/2;yi=y(i,i);K(1,2)=h*(

25、eval(fm);%计算 K2x=x;y(1,1)=y1+K(1,2)/2-K(1,1)/2;yi=y(i,i);K(1,3)=h*(eval(fm);%计算 K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;yi=y(i,i);K(1,4)=h*(eval(fm);%计算 K4y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4)/6; y(1,1)=y11(k1);x=a+(k1-1)*h;endy11else%以下计算过程为高阶初值问题mm=(b-a)/h;%共要求解mm个数据点for k2=1:m%读取初值条件fpri

26、ntf(请输入 %d 阶导数的初值 f(%d)(a)=n,(k2-1),(k2-1);y(k2,1)=i nput(=);endfor k2=1:my22(1,k2)=y(k2,1);%先把初值保存在矩阵y22(m, n)中,m表示第几个所求点，n表示第n阶初值endx=a;for k4=2:(mm+1)%求解mm个数据点的循环for k=1:(m-1)%计算K1，包括每一阶的 K1K(k, 1)=h*y(k+1,1);%y(k+1,1 )中 k+1 表示第 k+1阶，1表示第一个点；K(k, 1)中k表示阶数，1表示K1endK(m,1)=h*(eval(fm);x=x+h/2;%求解K1之

27、前，先重新对x和y赋值for k3=1:my(k3,1)=y(k3,1)+K(k3,1)/2;endfor k=1:(m-1)%计算 K2K(k, 2)=h*y(k+1,1);endK(m,2)=h*(eval(fm);x=x;for k3=1:my(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;endfor k=1:(m-1)%计算 K3K(k, 3)=h*y(k+1,1);endK(m,3)=h*(eval(fm);x=x+h/2;for k3=1:my(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2;endfor k=1:(m-1)%计算 K4K(k, 4)=h*y(k+1,1);endK(m,4)=h*(eval(fm);for k5=1:my22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4)/6;endfor k6=1:m%对 y(1,1), y(2,1)等重新赋值y(k6,1)=y22(k4,k6);endx=a+(k4-1)*h;endy22(:,1)end