数值代数与计算方法上机作业
（完整版）数值计算方法上机实习题答案
（完整版）数值计算⽅法上机实习题答案1．设?+=105dx xx I nn ，（1）由递推公式nI I n n 151+-=-，从0I 的⼏个近似值出发，计算20I ；解：易得：0I =ln6-ln5=0.1823, 程序为：I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为：20I = -3.0666e+010 （2）粗糙估计20I ，⽤nI I n n 515111+-=--，计算0I ；因为 0095.056 0079.01020201020≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为：I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。
⾸先分析两种递推式的误差；设第⼀递推式中开始时的误差为000I I E '-=，递推过程的舍⼊误差不计。
并记nn n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。
因为=20E 20020)5(I E >>-，所此递推式不可靠。
⽽在第⼆种递推式中n n E E E )51(5110-==-=Λ，误差在缩⼩，所以此递推式是可靠的。
出现以上运⾏结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2．求⽅程0210=-+x e x的近似根，要求41105-+?<-k k x x ，并⽐较计算量。
（1）在[0，1]上⽤⼆分法；程序：a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果：c =0.0903（2）取初值00=x ，并⽤迭代1021x k e x -=+；程序：x=0; a=1;while abs(x-a)>5*1e-4 a=x;x=(2-exp(x))/10; end x结果：x =0.0905（3）加速迭代的结果；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果：x =0.0995（4）取初值00=x ，并⽤⽜顿迭代法；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x结果： x =0.0905（5）分析绝对误差。
《计算方法》上机实验试题
《计算方法》上机实验试题
1. (25分)计算积分
dx x x I n n ⎰
+=105， n=0,1,2,…,20 若用下列两种算法 (A) n I I n n 151+
-=- (B) ⎪⎭
⎫ ⎝⎛-=-n n I n I 1511 试依据积分I n 的性质及数值结果说明何种算法更合理。
2. (25分)求解方程f(x)=0有如下牛顿迭代公式
)()(111---'-
=n n n n x f x f x x ， n ≥1，x 0给定 (1) 编制上述算法的通用程序，并以ε<--1n n x x （ε为预定精度）作为终止迭
代的准则。
(2)
利用上述算法求解方程 0cos :)(=-=x x x f
这里给定x 0=π/4，且预定精度ε=10-10。
3. (25分) 已知)3()(x x e x e x f -=，
(1)
利用插值节点x 0=1.00，x 1=1.02，x 2=1.04，x 3=1.06，构造三次Lagrange 插值公式，由此计算f(1.03)的近似值，并给出其实际误差； (2) 利用插值节点x 0=1，x 1=1.05构造三次Hermite 插值公式，由此计算f(1.03)的近
似值，并给出其实际误差。
4. (25分) 利用Romberg 算法计算积分
⎰
+4802cos 1dx x
精确到10
－4
总体要求：打印各题的程序代码及数值结果。
数值代数上机实验报告
数值代数课程设计实验报告姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法 二、实验环境 MATLAB7.0实验一、平方根法与改进平方根法一、实验要求：用熟悉的计算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解下列方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序分别求解40、84、120阶方程组的解。
二、算法描述及实验步骤GAuss 程序如下：（1）求A 的三角分解：LU A =;（2）求解b y =L 得y ； （3）求解y x =U 得x ；列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =;2求解b y P L =得y ； 3求解y x =U 得x ；三、调试过程及实验结果：%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7;>> for i=2:39b1(i)=15;end>> b1(40)=14;>> for i=2:83b2(i)=15;end>> b2(40)=14;>> for i=2:119b1(i)=15;end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1')>> x12=GAuss Zhu(A1,b1')>> x21=GAuss(A2,b2')>> x22=GAuss Zhu(A3,b3')>> x31=GAuss(A3,b3')>> x32=GAuss Zhu(A3,b3')运行结果：（n=40）GAuss消元法的解即为x11 =1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序：function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改进平方根法一、实验要求：用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、算法描述及实验步骤：平方根法函数程序如下：1、求A 的Cholesky 分解：L L A T=;2、求解b y =L 得y ；3、求解y x =TL 得x ； 改进平方根法函数程序如下：1、求A 的Cholesky 分解：T=LDL A ; 2、求解b y =L 得y ； 3、求解y x =TDL 得x ；三、调试过程及实验结果：clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11;>> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果：平方根法的解即为 x1 =1.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000改进平方根法解得的解即为x2 =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00090.99081.09080.1010四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend%------------------前代法求解Ly=b----------------------------for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法求解L'x=y-----------------------------A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法-----------------------------A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改进平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 -0.75 -0.5 0 0.25 0.5 0.75iy 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125i二、算法描述及实验步骤：QR分解求解程序如下：1、求A 的QR 分解;2、计算b c 11T =Q ;3、求解上三角方程1c x =R 得x ；三、调试过程及实验结果：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图（略）.（3）编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小，只需利用极值的必要条件0J a ∂=∂，0J b ∂=∂，0J c∂=∂，得到关于,,a b c 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 分别对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===，，，输入下列程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =0.3081 0.8587 1.4018 f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =0.3081 0.8587 1.4018f =924/2999*x^2+10301/11996*x+4204/2999>>。
上机数值计算练习题及答案.docx
习题31、在MATLAB 中,先运行指令A=magic(3), B=[l,2,l;3,4,3;5,6,7], C=reshape(l:6,3,2)生成阵列A 珂，B3X2,C3X2 ,然后根据运行结果回答以下问题：(1)计算A*B,B*A,这两个乘积相同吗？(2)计算A\B, B/A,左除、右除结果相同吗？(3)计算B(:,[1,2]).*C和C.*B(:, [1,2]),这两个乘积相同吗?(4)计算A\A和AAA,这两个计算结果相同吗？(5)计算A\eye(3)和inv(A),这两个计算结果相同吗？(提示：根据对计算结果的目测回答问题)解答如下：运行指令：A=magic(3), B=[l,2,l;3,4,3;5,6,7], C=reshape(l:6,3,2)得到结果：8 1 63 5 74 9 2B =1 2 13 4 35 6 7C =1 42 53 6(1)计算A*B,并得到结果如下：» A*Bans =41 56 5353 68 6741 56 45计算B*A, 并得到结果如下：»B*Aans =18 20 2248 50 5286 98 86从以上计算结果可以得出结论：A*BJ (2)计算A\B ，并得到结果如下：» A\Bans =0.0333 0.1000 0.16110.5333 0.6000 0.74440.0333 0.1000 -0.1722计算B/A, 并得到结果如下：B/Aans =0.0056 0.0889 0.17220.1389 0.2222 0.30560.2333 0.7333 0.2333 与B*A这两个乘积不同。
从以上计算结杲可以得出结论：左除、右除结杲不同。
(3)计算B(:,[1,2]).弋，并得到结果如下：A =» B(:,[1,2]).*C ans =1 8 6 20 15 36计算C.*B(:, [1,2]),并得到结果如下: » CFB(:, [1,2]) ans =1 6 20 15 36从以上计算结果可以得出结论:B(: J1,2]).*C 和C ・*B(:, [1,2])的两个乘积相同。
数值计算方法上机题目资料
2. 程序输入、输出用文件形式。 3. 编程语言要求用C, 编程环境TC或VC 4. 程序要求调试通过。 5. 每个方法要求给出一个具体的算例（可选
作业题）来验证。
五、上机报告要求
1.报告内容包括：
每种方法的算法原理及程序框图。 程序使用说明 具体算例及结果
上机调试体会及收获。
2.报告要手写。
六、上机报告及源程序提交时间
1.上机报告在考试当天提交。 2.源程序在考试前提交。
提交格式：文件夹（班级+姓名）
输入文件 程序文件夹 输出文件
不要拷贝其它文件！！！ 源程序
六、上机报告及源程序提交时间
源程序提交： 把以上文件压缩后，发送到以下邮箱：
haoyq@
七、考核方式
1.算法手算笔试（80%）+上机内容笔试 （10%）+上机报告（10%）
2.上机内容笔试可能形式：
编一段算法程序 给出一段算法程序，说明算法的名称。 程序填空 程序改错（包括算法和语法的错误）
数值计算方法上机练习
一、上机练习目的
复习和巩固数值计算方法的基本数学模型， 全面掌握运用计算机进行数值计算的具体 过程及相关问题。
利用计算机语言独立编写、调试数值计算 方法程序，培养学生利用计算机和所学理 论知识分析解决实际问题的能力。
二、上机练习任务
• 利用计算机基本C语言编写并调试一系列 数值方法计算通用程序，并能正确计算给 定题目，掌握调试技能。
• 掌握文件使用编程技能，如文件的各类操 作，数据格式设计、通用程序运行过程中 文件输入输出运行方式设计等。
• 写出上机练习报告。
三、数值计算方法上机题目
数值计算方法上机作业
《数值计算方法》上机作业学院：机械学院专业：机械设计及理论姓名：陈国保学号：2010412118日期：2010年12月18日第二章：插值法 2.1，问题：1.编制通用程序，对n+1个节点i x 及()ii y f x = (0,,)i n = （1）n 次拉格朗日插值计算公式Ln(x)； （2）n 次牛顿向前插值计算公式； （3）n 次牛顿向后插值计算公式； 2.已知()ln ,[,][1,2]f x x a b ==，取0.1,1,01,,10.i h x i hi ==+=用通用程序计算ln1.54及ln1.98的近似值。
流程图：2.2，源程序：主调函数：#include <stdio.h>#include <math.h>float Ln(float x1,int n,float x[80],float y[80]);float Nnf(float x1,int n,float x[80],float y[80]);float Nnr(float x1,int n,float x[80],float y[80]);main(){float x[80],y[80],x1,h,f_Ln,f_Nnf,f_Nnr;int n,i;//输入插值点printf("Please enter the interpolating point x:");scanf("%f",&x1);//输入拟合阶数printf("Enter the order n:");scanf("%d",&n);//输入已知数表的第一点x坐标printf("Enter the start point:");scanf("%f",&x[0]);//输入插值步长printf("Enter the step size h:");scanf("%f",&h);//计算已知数表的第一点y坐标y[0]=log(x[0]);//计算其余点的坐标for(i=1;i<=n;i++){x[i]=x[0]+i*h;y[i]=log(x[i]);}//拉格朗日法插值求解f_Ln=Ln(x1,n,x,y);//牛顿向前插值法求解f_Nnf=Nnf(x1,n,x,y);//牛顿向后插值法求解f_Nnr=Nnr(x1,n,x,y);//输出结果printf("The appoximate value caculated using Lagrange interpolation mathod is:\n%f\n",f_Ln);printf("The appoximate value caculated using Netown forward interpolation mathod is:\n%f\n",f_Nnf);printf("The appoximate value caculated using Netown backward interpolation mathod is:\n%f\n",f_Nnr);}拉格朗日插值函数：#include <stdio.h>#include <math.h>float Ln(float x1,int n,float x[80],float y[80]){int i,j;float sum=0;float nom=1,denom=1; //nom为基函数的分子，denom为基函数的分母 for(i=0;i<=n;i++){nom=1;denom=1;for(j=0;j<=n;j++)if(i==j) {nom*=1;denom*=1;} //当i=j时，分子分母不乘任何项 else {nom=nom*(x1-x[j]);denom=denom*(x[i]-x[j]);}//否则，分子自乘x-xj，分母自乘xi-xjsum=sum+nom/denom*y[i];//sum自加(x-xj)/(xi-xj)*yj}return sum;}牛顿向前插值函数：#include <stdio.h>#include <math.h>float Nnf(float x1,int n,float x[80],float y[80]){int i,j;float sum=0,t,h;float dif[80][80];float nom=1,denom=1;//nom为分子t(t-1)...(t-n+1)//denom为分母(n+1)!//计算步长hh=(x[n]-x[0])/n;//计算tt=(x1-x[0])/h;//构造向前差分表for(i=0;i<=n;i++)dif[0][i]=y[i];for(i=1;i<=n;i++)for(j=0;j<=n-i;j++)dif[i][j]=dif[i-1][j+1]-dif[i-1][j];sum=y[0];//计算插值公式for(i=1;i<=n;i++){nom=1;denom=1;for(j=1;j<=i;j++){nom*=(t-j+1);denom*=j;}sum+=dif[i][0]*nom/denom;}return sum;}牛顿向后插值函数：#include <stdio.h>#include <math.h>float Nnr(float x1,int n,float x[80],float y[80]) {int i,j;float sum=0,t,h;float dif[80][80];float nom=1,denom=1;//nom为分子t(t+1)...(t+n-1)//denom为分母n!//计算步长h=(x[n]-x[0])/n;//计算tt=(x1-x[n])/h;//构造向后差分表for(i=0;i<=n;i++)dif[0][i]=y[i];for(i=1;i<=n;i++)for(j=n;j>=i;j--)dif[i][j]=dif[i-1][j]-dif[i-1][j-1];sum=y[n];//计算插值公式for(i=1;i<=n;i++){nom=1;denom=1;for(j=1;j<=i;j++){nom*=(t+j-1);denom*=j;}sum+=dif[i][n]*nom/denom;}return sum;}计算结果：Please enter the interpolating point x:1.54Enter the order n:10Enter the start point:1Enter the step size h:0.1The appoximate value caculated using Lagrange interpolation mathod is:0.431782The appoximate value caculated using Netown forward interpolation mathod is: 0.431782The appoximate value caculated using Netown backward interpolation mathod is: 0.431782Press any key to continuePlease enter the interpolating point x:1.98Enter the order n:10Enter the start point:1Enter the step size h:0.1The appoximate value caculated using Lagrange interpolation mathod is:0.683097The appoximate value caculated using Netown forward interpolation mathod is: 0.683097The appoximate value caculated using Netown backward interpolation mathod is: 0.6830972.3计算结果分析：由以上计算结果可以看出采用10个点插值所得出的结果比较理想，误差很小，并且拉格朗日法，牛顿向前插值法和牛顿向后插值法计算的结果差不多。
数值上机作业1
数值上机作业1：Guass 顺序消去法信科08-2 0811620235 沈春燕一、小例题解方程组1231471258136101x x x ⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪= ⎪⎪ ⎪⎪⎪ ⎪⎝⎭⎝⎭⎝⎭程序：a=[1,4,7;2,5,8;3,6,10];b=[1,1,1]';n=3;for k=1:nfor i=k+1:nm(i,k)=a(i,k)/a(k,k); for j=k+1:na(i,j)=a(i,j)-m(i,k)*a(k,j); endb(i)=b(i)-m(i,k)*b(k); end endx(n)=b(n)/a(n,n); for i=(n-1):(-1):1 sum=0;for j=(i+1):nsum=sum+a(i,j)*x(j); endx(i)=(b(i)-sum)/a(i,i); end x运行结果：x =-0.3333 0.3333 0解的分析：通过手算可得出此题的精确解为11033T⎛⎫-⎪⎝⎭，通过上面的程序求解得出()0.33330.33330Tx =- ，此解基本符合精确解。
同时我还通过求出的解来反求向量b 时，我发现()111Tb =，与题中给出的b 相同，证明在求解阶数较小的时候，Guass 消去法是比较好的。
二、书本P39上机习题1：求解方程组12383846178611586115*861158614xxxxx⎛⎫⎛⎫⎛⎫⎪⎪ ⎪⎪⎪ ⎪⎪⎪ ⎪=⎪⎪ ⎪⎪⎪ ⎪⎪⎪ ⎪⎪⎪ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭⎝⎭程序：n=84;for l=2:n-1b(l)=15;endb(1)=7;b(n)=14;for i=1:na(i,i)=6;endfor j=1:n-1a(j,j+1)=1;endfor k=2:na(k,k-1)=8;endfor k=1:n-1for i=k+1:nm(i,k)=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m(i,k)*a(k,j); endb(i)=b(i)-m(i,k)*b(k);endendx(n)=b(n)/a(n,n);for i=(n-1):(-1):1sum=0;for j=(i+1):nsum=sum+a(i,j)*x(j);endx(i)=(b(i)-sum)/a(i,i);endx运行结果：x =1.0e+008 *Columns 1 through 60.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000Columns 7 through 120.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000Columns 13 through 180.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000Columns 19 through 240.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000Columns 25 through 300.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000000 0.000000010000001Columns 31 through 360.000000009999999 0.000000010000002 0.000000009999995 0.000000010000010 0.000000009999981 0.000000010000038Columns 37 through 420.000000009999924 0.000000010000153 0.000000009999695 0.000000010000610 0.000000009998779 0.000000010002441Columns 43 through 480.000000009995117 0.000000010009765 0.000000009980470 0.000000010039060 0.000000009921880 0.000000010156240Columns 49 through 540.000000009687519 0.000000010624962 0.000000008750076 0.000000012499847 0.000000005000305 0.000000019999390Columns 55 through 60-0.000000009998779 0.000000049997559 -0.000000069995117 0.000000169990233 -0.000000309980464 0.000000649960918Columns 61 through 66-0.000001269921799 0.000002569843445 -0.000005109686279 0.000010249370117 -0.000020468730470 0.000040967421880Columns 67 through 72-0.000081904687519 0.000163838750076 -0.000327645000305 0.000655310001221 -0.001310550004883 0.002620970019531Columns 73 through 78-0.005241270078125 0.010480010312500 -0.020949751250000 0.041858575000000 -0.083553290000000 0.166451289999999Columns 79 through 84-0.330281269999999 0.650077449999998 -1.258214389999996 2.348666889999992-4.026286069999986 5.368381449999981解的分析：当n 从1到40时，求出的解全为1，没有任何的误差；当n=41时，求出的解出现了较为微小的误差，先是只有41x 以后出现很小的误差，后来随着n 的增加，逐渐的35x 至41x 之间也逐渐出现了误差并且误差随着n 的增加而逐渐增加。
数值计算方法上机实验
用追赶法解三角方程组Ax=d，A= ， d=
实验过程：
1.利用C语言编程求解，程序如下:
#include<stdio.h>
#include<math.h>
void main()
{int i;
double a[5]={0,-1,-1,-1,-1},b[5]={2,1,1,1,1},c[4]={2,2,2,2},f[5]={6,7,9,11,1};
{double x=-0.99,y;
int n=0;
printf("%d ,%lf\n",n,x);
do
{y=x;
x=x-(x*x*x/3-x)/(x*x-1);
n++;
printf("n= %d ,x= %lf\n",n,x);
}while(fabs(y-x)>1e-5);
}
1.2.命令窗口结果截屏如下：
}
2.2.命令窗口结果截屏如下：
实验结果分析：
从运算截图中我们可以知道，一方面运用牛顿迭代法和牛顿下山迭代法的程序运算的结果是一致的，若出现小小的差异的话，可能是由于设置的变量类型不同导致，另一方面运用牛顿下山迭代法比运用一般牛顿迭代法计算的步骤要少一些，这证明对于同一个题，牛顿下山迭代法的优越性比一般的牛顿迭代要高。
实验内容：
用列主元法解线性方程组
=
实验过程：
1.利用C语言编程求解，程序如下：
#include<stdio.h>
#include<math.h>
void main()
{double t,m1,m2,m3,x1,x2,x3,x4,a[4][5]={1.1348,3.8326,1.1651,3.4017,9.5342,0.5301,1.7875,2.5330,1.5435,6.3941,3.4129,4.9317,8.7643,1.3142,18.4231,1.2371,4.9998,10.6721,0.0147,16.9237};
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值代数与计算方法上机作业
作业一：Matlab的基本操作
P31
1.根据习题12和习题13构造算法和MATLAB程序，以便精确计算所有情况下的二次方程的根，包括b?b2?4ac的情况。
2.参照例1.25，对下列3个差分方程计算出前10个数值近似值。
在每种情况下引入一个笑
?1?得出是误差。
如果没有初始误差，则没个差分方程将生成序列?n? 。
构造类似表1.4、
2??n?1表1.5以及图1.8至图1.10的输出。
?1rn?1，其中n=1,2,… 23(b) p0?1,p1?0.497,pn?pn?2, 其中n=2,3,…
25(c) q0?1,q1?0.497,qn?qn?1?qn?2, 其中n=2,3,…
2(a) r0?0.994；rn?
作业二：非线性方程f(x)?0的解法
P40
1. 使用程序
2.1求解下面每个函数的不动点（尽可能多）近似值，答案精确到小数点后12为。
同时，构造每个函数和直线y=x来显示所有不动点。
（a）（b）（c）（d）
g(x)?x5?3x3?2x2?2
g(x)?cos(sin(x)) g(x)?x2?in(x?0.15) g(x)?xx?cos(x)
P49
3. 修改程序2.2和程序2.3，使得输出分别类似于表2.1和表2.2的矩阵（即矩阵的第一行应当为[0
a0 c0 b0 f(c0) ] ）。
P69 4，用习题11中的立方根算法修改程序2.5，并用其近似下列每个立方根到小数点后10位。
（a） p0?2，求7的近似值。
（b） p0?6，求200的近似值。
（c）
p0??2，求(?7)的近似值。
作业三：线性方程组AX?B的求解方法
131313P93 1. P97 2.
P109 2.
P120 1.
P130 4.
作业四：插值与多项式逼近
P154
Matlab的矩阵特性使其能够快速计算一个函数在其多个点处的值。
例如，如果X=[-1 0 1],则sin(X)将得到[sin(-1),sin(0),sin(1)]。
类似地，如果X=-1:0.1:1，则Y=sin(X)将得到与X同样维数的矩阵Y,其值为正弦函数的值。
通过定义矩阵D=[X’，Y’],可将这两个行矩阵输出为表的形式。
注意：矩阵X和Y必须有相同的长度。
1. （a）用plot命令，在同一幅图正绘制区间-1?x?1上的sin(x),习题1中计算出的P5(x)，
P7(x)和P9(x)。
（b）创建一个表，他的各列分别由区间[-1,1]上的10个等距点x处的sin(x)，
P5(x)，P7(x)和P9(x)值构成。
P160
1. 用Matlab实现算法4.1,多项式P(x)?aNxN?aN?1xN?1???a2x2?a1x?a0的系数以
1?N矩阵P?[aN,aN?1,?,a2,a1,a0]的形式输出。
P171
2. 下表给出了11月8日美国洛杉矶的一个郊区在5小时内的测量温度。
（a）利用程序4.1，对表中的数据构造一个拉个人朗日插值多项式。
（b）利用算法4.1（iii）,
估计在这5小时内的平均温度。
（c）在同一坐标系中画出表中的数据和由(a)得到的多项式。
讨论用(a)中的多项式
计算平均温度可能产生的误差。
时间（下午）华氏度
1 66
2 66
3 66
4 64
5 63
6 63
P178
1,用程序4.2重新计算4.3.5节中的第2题
1.胡克（Hooke）定律指出F?kx,其中F是拉伸弹簧的拉力（单位为盎司），x是拉伸的长度（单位为英寸）。
根据下列试验数据，实用程序5.1求解拉伸长量k的近似值。
（a）（b）见书上
作业五：数值微分
P260
1, 用程序6.1求解下列函数在x处的导数近似值，精度为小数点后13位。
注：有必要改写程序中的max1的值和h的初始值。
（a） f(x)?60x45?32x33?233x5?47x2?77; x?1/3
（b） f(x)?tan(cos(5?sin(x)1?5));x? 231?x（c） f(x)?sin(cos());
1xx?12
（d） f(x)?sin(x?7x?6x?8); （d） f(x)?x;
P270
xx32x?1?5 2x?0.0001
1. 修改程序6.3，使`得可用它计算P?(xM),M?1，2，?，N?1
作业六：数值积分
P290
1, (a) 对习题1中的每个积分，计算M和步长h,使得用组合梯形公式计算得到精确到小数
点后9位的结果。
用程序7.1计算每个积分。
（b）对习题1中的每个积分，计算M和步长h,使得用组合辛普森公式计算得到精确到小数点后9位的结果。
用程序7.2计算每个积分。
P301
1, 利用程序7.4求习题1中的积分，精确到小数点后11位。
P307
1, 用程序7.6求以下定积分的近似值，其实容错?0?0.00001。
（a）（b）
?30sin(2x)dx 51?x??301sin(4x)e?2xdx
（c）
0.04
高斯赛德尔迭代算法代码：function X = gseid(A, B, P , delta, max1) N = length(B); for k=1:max1 for j =1:N if j==1 X(1)=(B(1)-
A(1,2:N)*P(2:N))/A(1,1); elseif j==N X(N)=(B(N)-A(N, 1:N-1)*(X(1:N-1))')/A(N, N); else X(j)=(B(j)-A(j, 1:j-1)*X(1:j-1)'-A(j,
j+1:N)*P(j+1:N))/A(j,j); end end err=abs(norm(X'-
P));relerr=err/(norm(X)+eps);P=X'; if(err> B = matrix2(50)>> A =
matrix(50)>> X = gseid(A, B, B, 0.0001, 50)X = 0.4638 0.5373 0.5090 0.4982 0.4989 0.5000 0.5001 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5001 0.5000 0.4989 0.4982 0.5090 0.53730.4638
114569384
感谢您的阅读，祝您生活愉快。