(完整word版)计算方法上机题目(2)

1、目录1. 计算方法 A 上机作业 1.上机练习目的 1.上机练习任务 1.计算方法 A 上机题目 1.程序设计要求 1.上机报告要求 1.2. QR 分解法求解线性方程组 2.计算原理 2.程序框图 7.计算实习 8.Matlab 代码 8.3. 共轭梯度法求解线性方程组 1.0计算原理 1.0.程序框图 1.1.计算实习 1.2.Matlab 代码 1.2.4. 三次样条插值 1.4.计算原理 1.4.程序框图 1.6.计算实习 1.7.Matlab 代码 1.7.5. 四阶龙格库塔法求解常微分方程的初值问题 2. 1计算原理 2.1.程序框图 2.2.计算实习 2.3.Matlab 代码

2、2.3.1. 计算方法 A 上机作业上机练习目的复习和巩固数值计算方法的基本数学模型，全面掌握运用计算 机进行数值计算的具体过程及相关问题。利用计算机语言独立编写、调试数值计算方法程序，培养学生 利用计算机和所学理论知识分析解决实际问题的能力。上机练习任务? 利用计算机语言编写并调试一系列数值方法计算通用程序，并 能正确计算给定题目，掌握调试技能。? 掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。? 写出上机练习报告。计算方法 A 上机题目1. QR 分解方法求解线性方程组。 （第二章）2. 共轭梯度法求解线性方程组。 （第三章）3. 三次

3、样条插值（第四章）4. 四阶龙格库塔法求解常微分方程的初值问题程序设计要求1. 程序要求是通用的，在程序设计时要充分考虑哪些变量应该可 变的。2. 程序要求调试通过。上机报告要求报告内容包括：每种方法的算法原理及程序框图。程序使用说明。算例计算结果。2. QR分解法求解线性方程组计算原理x vPx v F2当x Rn是任意给定的非零向量，v Rn是任意给定的单位向量，则存在初等反射阵H I 2uuT ,使得Hxv,其中 为常数，当取单位向量u时，由u确定的矩阵H必定满足Hxv ,所以在计算过程中取u的值为上述值。设A是一个m n m n阶矩阵且它的列向量线性无关，则利用豪斯霍尔德变换可以把A逐

4、步化为上梯形矩阵，设ana21AMam1a12a22Mam2a1na2nMamna1, a2, L , an7具体变换过程如下:1,2,L ,n是m维单位坐标向量。为把矩阵A的第一列a10化为1,0,L ,0 T，取0,v1,0丄T,0Rm （或取v e1），根据上式可得，取0a1Pa11e pu1-0-其中0Pai1e0aiie2uuT4P 1 P201 a11T1 1令H11 m2u1u1* 1 m2 1 111T，用H1左乘A得,A1H1A00 0H 1a1 , H 1a2 ,L1042 1 a11卫彳1丄t 6昂0hV绘*B 0HAam2A cm34盘/其中曙=局硏 =（爲- 口=曙

5、- jSijWi = a畀- 0丄,（af）- 66）f曙-旳（后-6）4T 一 ftjOa?由此可知=增_為衆_的）.j = 2A侃4P =谓 一厲皿器，i = 2,驚,m;j = 2,3,.坐标向量.这时，取乃-一桦1（&）|詁“应.第2步记铲亠（翻严朗，观）丁 -（翻应严”显然蹭丰0* av = f）（或取v = 一屮、其中411为切=（0詁叽）T中的w-1维耳其中II讪b =忆0 -叱Ph = *2喝-卩埸）,2蚀u成wqwJ卩如临血伽-國?）其中构造令肓2 = /m_i - 2u2U2 = An-1 - 計沖孑,Z)用血左乘且得4卩）=地&=（6 6&；）引+帀切严件，血）d盘丿次变

6、换，即存在初轩反射阵g.a；气 *HI-6第上步记。）=（讹堞，总HP出Q）T=（4：）学,，a匕?,址 7）T,琶i =（硝7,监呂，话】）T, j = &上+1，.显然“旷】）* 0,取丄=av =占“）（或取乜=-皆7）,其中才7为 s =中的刃一上+ 1维单位坐标向量.这时取Vr1个.Wk憚7_叶占7h l|wfe|25其中6 = -sgn网：7)|吐i|b, wfc = 7 - 叶ef而II讥lb = |进7 - trfc4k_l)lb =丿2(來-%&：-】)，T2wkuwkwWfcW?2吸 _ WF 一叶仇吐：7)=匚厂，其中otjc =丑（乐一心：“）.构造片At = An亠*

7、+1 f= Im-k+l 3丿七 j令用Hk左乘Ad)得A( = fffc4(*_1) =HA(k-i65,a诙i +巾也,，嚟勺+ 4斑,1注1卫0j=i1 6%仪al5fc+ll.i+2g 嚟7 0(3)a2k*fi2A+l*(2)a2.k+2-* :Ji-nVi点I)叫花41rt(*)%住+2(k)k+】我+1(fc)4+12 *t*如其中(執疋护=月站1J = (/m JH 1 -叫1T (fc = 1叫応)a)7由此即得Oy =硝-。-%泌：7如-讥+丈 盘旷心 mi“Xi-*+l)继续上述变换过程，直至9二minm - 1加步，懾后得到A =H“= R*由于每个 乐 都是正交矩阵，

8、所以HH、也是正交矩阵上式可以写为HA 工 &(246)当m = n时，R是n阶上三角矩阵-这时.式(246)可以写为A=QR,(247)其中Q = Wr为n阶正交矩阵，即把n阶矩阵A分解为一个n阶正交矩阵与一 个可逆的n阶上三角矩阵R的乘积程序框图，/输入系数矩阵 A/ 右端列向量bx8计算实习用豪斯霍尔德变换求方程组Axb其中54756753941287886537810987756A57911975,b 5368891089587877810105756759101052Matlab 代码%使用说明：%需自己输入矩阵 A 及 b 的值%变量Q,R分别为进行QR分解后的结果clearclc

9、format longload(A 矩阵 .mat)load(b 矩阵 .mat)%调用函数 qrhs 进行 QR 分解Q,R=qrhs(A);,n=size(R);fprintf( 您输入的矩阵阶数 );ny=Q*b;%回代过程x(n,1)=y(n,1)/R(n,n);for i=n -1:-1:1s=y(i,1);for j=i+1:ns=s-R(i,j)*x(j,1);end x(i,1)=s/R(i,i);endx被调用函数 qrhs functionQ,R=qrhs(A) format long ,n=size(A);Q=eye(n);for j=1:n -1B=norm(A(j:n

10、,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=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)，由

11、迭代格式kd产生的迭代序列在无舍入误差的假定下，最多经过n次迭代就能求得f x的最小点，也就是方程组Ax b的解。共轭梯度法中关键的两点是，确定迭代格式中的搜索方向和最佳步长。 搜索方向d k，与前一次的搜索方向关于d k1关于矩阵A共轭，即dk,Adk1 0，然后从点x，沿d k方向求得f x的极小值点x k 128minf由此解得最佳步长k和参数k的表达式为k T . k rk _ddk T Adk 1 T k r Ad,k T , k d Ad共轭梯度法的计算公式为:d0kx(k 1) rbrkTdk xAx0k.kAdA kk d“ k 1Ax1 TkAddk1,k T kd Ad(k

12、 1)krkd程序框图计算实习用共轭梯度法求解线性方程组Ax b ，其中2111210AO O O ,bM1 2 10121矩阵 A 的阶数 n 分别取为 100,200,400，指出计算结果是否可靠 Matlab 代码%使用说明：共轭梯度法求解 Ax=b %命令行中输入矩阵 A 及 b %然后调用函数 getd 进行计算%变量含义：n方程阶数，x0初始向量%计算精度，r残向量clearclcformat 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

13、;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 和 bx0=zeros(n,1);%生成初始向量e=input(请输入计算精度e=);%输入计算精度x=getd(A,b,x0,e);%调用函数 getd 进行计算if n=100%对 x 元素进行重新排列x=reshape(x,10,10)elseif n=200 x=reshape(x,10,20) else x=reshape(x,20,20) end 被调用函数 getdfunction x=getd(A,b,xO

14、,e)% 矩阵 A,b，初始向量 x0，精度 e n=size(A,1);% 获取矩阵 A 的阶数x=xO;%初始向量 r=b-A*x;% 残向量 d=r;%搜索向量for k=O:(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=n orm(r2F2/norm(rF2;d=r2+q*d;r=r2;end4. 三次样条插值计算原理设在区间4方上绪定】个节点厲ax xi - xb) ,tE节点斗处的函数 值为yi -Xxr) (/若函数S(x)満足以下三个条件】(1)在每*子区间EUo丄期匕 乂

15、4是三次多项式1 Sfo) Vf (i -G)在b上.S(Q的二阶导数S(r)连绽.世垛SM为甫数胁-JW在区闾e b上的三次祥条插值函数由定义可知SC护有4两个持定鑫数.根据条件(3)可得如下3併3心方程.s：00 $；(%),(5)血)畑，再加上条件(2)的卄1亍方程.共有 X 个方程*因此还需嬰增如两个亲件才可以确定 出如个特定参数.所増J0的条件祢为边界条件或划点条件常用的三种边界条件为:已知厂(“厂0”已/T(a)-r()i己知冗0是以卩g为周期的周KSK.ft节点曲业的二阶导数值为Mt(i = 0J/-A M,为待定魯数眾据公式0)可以毎到区间&】珂上的表达式为要得到5(的表达式*

16、就要确定M的值.越于奉同的边界条件，求馨M的力理生 不同弟一类*痢二类和第三类边界条件的三弯矩方程组分别为公式(7X92 : 码2右* VV=41 川 Af d2(?)J 2匚1P-)3 )21冏2 *f4 f ,a申2中1 2比丿虬丿阳M,9丿W4Vdy&2K 2(9)其中: I * =1 / = 11 -1.aj公式中的4 = 6/和斥坷n.=6”/J公式(9)中的町严L心“如 耳+州九钠I丿三次特条插值逓数的关權是3E解式(7). (8)(9示的方程齟.将得到的 制代入 公武3)中可帶判相应的插值程序框图计算实习给定函数f(x) 1 2( 1 x 1).取等距节点，构造三次样条插值S10

17、(x) 1 15xMatlab代码%使用说明：该程序解决的是三次样条插值中，第1,2类边界条件的问题%各变量含义：a,b插值区间左右端点%n插值节点数目%p,q左右端点导数值%A，M，d用于求解AM=d中，矩阵M的值%C存放各区间内插值函数的系数矩阵%zglu利用追赶法进行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=zero

18、s( 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%选择边界条件进行计算，并输入区间左右端点的导数值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)=

19、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,:)

20、-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=d M=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,:)A2*M(k -1,:)+1/h*(y(k,:) -1/6*hA

21、2*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,:)-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 l=1:n-1disp(num2str(C(

22、l,1),xA3+,num2str(C(l,2),xA2+,num2str(C(l,3),x+,num2str(C(l,4),num2str(x(l,:), xtwW = 14r-* tVm(a) = y1程序框图计算实习用标准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)=yclearc

23、lcformat longm=i nput(请输入常微分方程的阶数m=);a=input(请输入x下限a=);b=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)y1=y(1,1);%计算 K1K(1,

24、1)=h*(eval(fm);x=x+h/2;y(1,1)=y1+K(1,1)/2;y1=y(1,1);K(1,2)=h*(eval(fm);%计算 K2x=x;y(1,1)=y1+K(1,2)/2 -K(1,1)/2;y1=y(1,1);K(1,3)=h*(eval(fm);%计算 K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;y1=y(1,1);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%读取初值条件fprintf(请输