计算方法上机实习题

合集下载

计算方法与实习答案

计算方法与实习答案【篇一：《基础会计学习指导、习题与实训》答案】名词解释1．会计：是以货币为主要计量单位，以凭证为依据，运用专门的技术方法，对一定主体的经济活动进行连续、系统、全面的核算与监督，以提高经济效益为目标，向有关方面提供会计信息的一种经济管理活动。

2．会计职能：是指会计在经济管理中所具有的功能，即会计在经济管理中能发挥什么作用。

3．会计核算职能：是指以货币为主要计量单位，对企事业单位一定时期的经济活动进行真实、连续、系统、完整的记录、计量和报告。

4．会计监督职能：是指依据监督标准，利用会计核算所提供的会计信息对各单位的经济活动全过程的合法性、合理性和有效性进行的指导、控制和检查。

5．会计对象：是指会计所要核算和监督的内容，即会计工作的内容。

6．会计要素：是对会计对象按经济特性所做的基本分类，是会计对象的具体内容。

7．资产：是指企业过去的交易或者事项形成的、由企业拥有或者控制的、预期会给企业带来未来经济利益的资源。

8．负债：是指企业过去的交易或者事项形成的、预期会导致经济利益流出企业的现时义务。

9．所有者权益：是指企业资产扣除负债后由所有者享有的剩余权益，包括实收资本、资本公积、盈余公积和未分配利润。

10．收入：是指企业在日常活动中形成的、会导致所有者权益增加的、与所有者投入资本无关的经济利益的总流入，包括销售商品收入、劳务收入、利息收入等。

11．费用：是指企业在日常活动中发生的、会导致所有者权益减少的、与向所有者分配利润无关的经济利益的总流出。

12．利润：是指企业在一定会计期间的经营成果，包括收入减去费用后的净额、直接计入当期利润的利得和损失等。

13．会计方法：是为实现会计核算、进行会计管理和完成会计任务所采用的手段。

14．会计核算方法：是对单位已经发生的经济活动进行连续、系统、全面的核算所采用的方法，包括设置账户、复式记账、审核和填制会计凭证、登记账簿、成本计算、财产清查和编制财务会计报告。

Access学习辅导与上机实习—习题(第3章查询)

Access学习辅导与上机实习—习题（第3章查询）■典型习题一、选择题1、是_________不会在简单查询的设计网络中出现。

A .字段B 条件C更新D排序2、下列叙述不正确的是。

A 删除查询主要用于删除符合条件的记录。

B 更新查询中可以使用计算功能。

C 追加查询时如果两个表的结构不一致，则不能进行。

D生成表查询生成新的表，该表是源表的一个子集。

3、不属于统计函数。

A MAXB COUNTC LASTD Y EAR二、填空题4、在查询表达式中，“/”表示，而“\”表示。

5、总计查询中，必须包含的两种字段是和。

6、在查询设计过程中，有多种方式可以视察查询结果，比如可以进入视图模式或单击按钮。

三、操作题7、对“教学”库建立下列查询：（1）用库中的三张表建立一个由学生姓名、课程名称、成绩组成的新数据表，新表取名为新表1。

（2）用学生表生成一个关于女生情况的新数据表，取名为女生情况表。

（3）用库中的三张表建立一个由学生姓名、成绩组成的数学成绩表，取名为数学成绩表。

（4）用库的三张表建立一个不及格学生的汇总表，表中的字段要求有学生姓名、课程名称、成绩。

（5）用库中的三张表建立一个由课程名称、各科平均分、总分、最高分、最低分组成的新数据表，取名为成绩统计表。

8、查询指定课程名的学生成绩，课程名只要求输入任意的几个字。

9、查询统计指定课程名称的总分和平均分。

■习题一、选择题1、利用对话框提示用户输入参数的查询过程称为_______。

A 选择查询B 参数查询C 操作查询D SQL查询2、以下叙述中，错误的是______。

A 查询是从数据库的表中行筛选出符合条件的记录，构成一个新的数据集合B 查询的种类有：选择查询、参数查询、交叉查询、操作查询和SQL查询C 创建复杂的查询不能使用查询向导。

D 可以使用函数、逻辑运算符、关系运算符创建复杂的查询3、Access共提供了种数据类型。

A 8B 9C 10D 114、用向导创建简单查询的步骤不包括。

北航数值分析计算实习题目二矩阵QR分解

数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。

已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。

说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。

2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)；矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。

2. 算法设计方案本题采用带双步位移的QR 分解方法。

为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。

在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。

同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===，其中k 为迭代次数。

Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4：如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。

东南大学计算方法与实习上机实验一

东南大学计算方法与实习实验报告学院：电子科学与工程学院学号：06A*****姓名：***指导老师：***实习题14、设S N=Σ(1)编制按从大到小的顺序计算S N的程序；(2)编制按从小到大的顺序计算S N的程序；(3)按两种顺序分别计算S1000，S10000，S30000，并指出有效位数。

解析：从大到小时，将S N分解成S N-1=S N-,在计算时根据想要得到的值取合适的最大的值作为首项；同理从小到大时，将S N=S N-1+ ，则取S2=1/3。

则所得式子即为该算法的通项公式。

(1)从大到小算法的C++程序如下：/*从大到小的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;const int max=34000; //根据第(3)问的问题，我选择了最大数为34000作为初值void main(){int num;char jus;double cor,sub;A: cout<<"请输入你想计算的值"<<'\t';cin>>num;double smax=1.0/2.0*(3.0/2.0-1.0/max-1.0/(max+1)),temps;double S[max];// cout<<"s["<<max<<"]="<<setprecision(20)<<smax<<'\n';for(int n=max;n>num;){temps=smax;S[n]=temps;n--;smax=smax-1.0/((n+1)*(n+1)-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/num-1.0/(num+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-smax); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<n<<"]="<<setprecision(20)<<smax<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(2)从小到大算法的C++程序如下：/*从小到大的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;void main(){int max;A: cout<<"请输入你想计算的数，注意不要小于2"<<'\t';cin>>max;double s2=1.0/3.0,temps,cor,sub;char jus;double S[100000];for(int j=2;j<max;){temps=s2;S[j]=temps;j++;s2+=1.0/(j*j-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/j-1.0/(j+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-s2); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<j<<"]="<<setprecision(20)<<s2<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(3)(注：因为程序中setprecision(20)表示输出数值小数位数20，则程序运行时所得到的有效数字在17位左右)ii.选择从小到大的顺序计算S1000、S10000、S30000的值需要计算的项S1000S10000S30000计算值0.74900049950049996 0.74966672220370571 0.74996666722220728实际精确值0.74900049950049952 0.74990000499950005 0.74996666722220373误差 4.4408920985006262*10-16 5.6621374255882984*10-15 3.5527136788005009*10-15有效数字17 17 17附上部分程序运行图：iii.实验分析通过C++程序进行计算验证采用从大到小或者从小到大的递推公式算法得到的数值基本稳定且误差不大。

实验六上机实习_2016

实验6 地下水动力学上机实习一、实习目的与要求（1）熟练掌握含水层试验软件AquiferTest Pro V2016的操作使用方法。

（2）基本掌握使用软件AquiferTest Pro V2016中的Theis 和Cooper-Jacob Time-Drawdown 分析方法进行抽水实验（Pumping Test ）求解含水层参数。

二、分析方法（1）Theis 方法（标准曲线配比法）在双对数坐标纸上，以W (u )为纵坐标，1/u 为横坐标作出的曲线通常称为泰斯曲线。

在双对数坐标纸上以t 或t /r 2为横坐标，s 为纵坐标点出观测数据，通过对比观测数据点与泰斯曲线来求解含水层参数。

（2）Cooper-Jacob Time-Drawdown （直线图解法）Cooper-Jacob (1946)方法是简化的泰斯方法，通过有效地增大时间值、减小井距（即减小u 值）来计算。

直线图解法将各个观测孔的降深随时间变化的数据标在半对数坐标纸上，以时间为横坐标，降深为纵坐标。

如果有足够多的数据就可以连成一条直线，含水层的导水系数和储水系数就可以按照以下两个方程计算出结果：s Q T ∆=π43.2，2025.2r Tt S = 三、软件简介AquiferTest V2016①整合了抽水试验Pumping Test 和微水试验Slug Test 数据分析技术，本次上机实习主要学习前者。

Aquifer Test 4.2具有友好的界面，快捷易用的优点，并且拥有计算各类含水层特性的功能（包括承压含水层、非承压含水层、越流含水层、裂隙含水层和井储等），对于同一份数据可以通过建立多种分析方法来做比较，并可生成专业的抽水试验报告，图件都可以图片形式导出。

试验所需数据可以通过键盘输入，也可以通过导入Microsoft Excel 文件或ASCII 格式文件。

对于抽水试验，Aquifer Test 2016提供以下14种解决方法：➢Time-Drawdown （时间-降深分析） ➢Time Drawdown-Discharge ➢Theis （承压含水层） ➢Theis 结合Jacob 修正(非承压含水层） ➢Neuman （非承压含水层） ➢Boulton （非承压含水层） ➢Hantush-Jacob （越流含水层，有存储的弱含水层） ➢Hantush （越流含水层，有存储的弱含水层） ➢Walton （双孔隙度，裂隙流） ① 目前的版本为AquiferTest Pro 2016版本，新版本可以通过https:///aquifertest/下载。

数值方法计算实习题

数值⽅法计算实习题数值⽅法计算实习题⼀、下表给出了飞⾏中鸭⼦的上部形状的节点数据，试⽤三次样条插值函数（⾃然边界条件）和20次Lagrange 插值多项式对数据进⾏插值。

⽤图⽰出给定的数据，以及()s x 和20()L x 。

12 12.6 13.0 13.3];>> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25]; %（1）三次样条插值法xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> title('试验⼀--三次样条插值图⽰')024********试验⼀--三次样条插值图⽰>> pp=spline(x,y)pp =form: 'pp'breaks: [1x21 double]coefs: [20x4 double]pieces: 20order: 4dim: 1>> pp.coefsans =0.7735 -0.9995 0.7760 1.3000 0.7735 -0.0714 0.3477 1.5000 -2.7894 1.3209 1.0974 1.8500 -0.4585 -0.3528 1.2910 2.10000.4489 -1.0405 0.5944 2.6000 0.1738 -0.5018 -0.0225 2.70000.0783 -0.0325 -0.5033 2.40001.3141 0.0850 -0.47712.1500 -1.5812 1.2676 -0.0713 2.0500 0.0431 -0.1555 0.2623 2.1000 -0.0047 -0.0261 0.0808 2.2500 -0.0245 -0.0401 0.0146 2.3000 0.0175 -0.1135 -0.1390 2.2500 -0.0128 -0.0505 -0.3358 1.9500 -0.0201 -0.1003 -0.5319 1.4000 1.2094 -0.1485 -0.7310 0.9000 -0.8279 0.9400 -0.4935 0.7000 0.0122 -0.0535 -0.1389 0.6000 -0.2960 -0.0316 -0.1900 0.5000 -0.2960 -0.3867 -0.3573 0.4000 所以所得⽅程为%（2）⽤拉格朗⽇法插值%定义Lagrange程序function f=Language(x,y,x0)syms t;if(length(x)==length(y))n=length(x);elsedisp('xoíyµêy2??àµè￡?');return;endf=0.0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end ;for (j=i+1:n)l=l*(t-x(j))/(x(i)-x(j)); end ; f=f+l; simplify(f); if (i==n)if (nargin==3) f=subs(f,'t',x0); elsef=collect(f); f=vpa(f,6); end end end>> Language(x,y) ans =52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2⼆、已知Wilson 矩阵1078775658610975910A=，且向量32233331b ??=，则⽅程组Ax b =有准确解[]1111Tx =。

实习题1

计算方法实验报告班级：学号：姓名：成绩：1 舍入误差及稳定性一、实验目的（1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；（2）通过上机计算，了解舍入误差所引起的数值不稳定性二、实验内容1、用两种不同的顺序计算∑n -2，分析其误差的变化n =1100002、已知连分数f =b 0+a 1，利用下面的算法计算f ：b 1+a 2/b 2+a 3/(...+a n /b n ) d n =b n , d i =b i +a i +1f =d 0 (i =n -1, n -2, ... , 0d i +1写一程序，读入n , b 0, b 1,..., b n , a 1,..., a n , 计算并打印f3、给出一个有效的算法和一个无效的算法计算积分x n1y n =⎰dx (n =0, 1, ... , 04x +114、设S N =∑j =2N 11⎛311⎫，已知其精确值为-- ⎪ 2j -12⎝2N N +1⎭（1）编制按从大到小的顺序计算S N 的程序（2）编制按从小到大的顺序计算S N 的程序（3）按两种顺序分别计算S 1000, S 10000, S 30000, 并指出有效位数三、实验步骤、程序设计、实验结果及分析1、用两种不同的顺序计算∑n ，分析其误差的变化n =110000-2（1）实验步骤：分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计：a. 顺序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=1;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sun[%d]=%-30f",n,sum); if(n>=10000)break; n++;}printf("sum[%d]=%f\n",n,sum);}b. 逆序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=10000;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sum[%d]=%-30f",n,sum);if(n<=1)break;n--;}printf("sum[%d]=%f\n",n,sum); }3）实验结果及分析：程序运行结果：a. 顺序计算b. 逆序计算结果分析：两种不同顺序计算结果是一样的，顺序计算误差从一开始就很小，而逆序计算误差最开始十分大，后来结果正确。

《计算方法》上机实验指导书刘轶中-8页word资料

理学院《计算方法》实验指导书适合专业：信息与计算科学数学与应用数学贵州大学二OO七年八月前言《计算机数值计算方法》包括很多常用的近似计算的处理手段和算法，是计算科学与技术专业的必修课程，为了加强学生对该门课程的理解，使学生更好地掌握书中的计算方法、编制程序的能力，学习计算方法课程必须重视实验环节，即独立编写出程序，独立上机调试程序，必须保证有足够的上机实验时间。

在多年教学实践基础上编写了《计算机数值计算方法》上机实习指南，目的是通过上机实践，使学生能对教学内容加深理解，同时培养学生动手的能力.本实习指南，可与《计算机数值计算方法》课本配套使用，但是又有独立性，它不具体依赖哪本教科书，主要的计算方法在本指南中都有，因此，凡学习计算方法课的学生都可以参考本指南进行上机实习。

上机结束后，按要求整理出实验报告。

实验报告的内容参阅《计算机数值计算方法》上机实验大纲。

目录第一章解线性方程组的直接法实验一 Gauss列主元素消去法实验二解三对角线性方程组的追赶法第二章插值法与最小二乘法实验三 lagrange插值法实验四分段插值法实验五曲线拟合的最小二乘法第三章数值积分实验六复合求积法实验七变步长法第四章常微分方程数值解法实验八 Euler 方法第五章解线性方程组和非线性方程的迭代法实验九 Jacobi 迭代法、Gauss-Seidel 迭代法实验十 Newton 迭代法实验一： Gauss 列主元素消去法实验学时：2实验类型：验证实验要求：必修一、实验目的用gauss 消去法求线性方程组AX=b. 其中一、实验内容二、实验条件PC 机，tc2.0，Internet 网。

三、实验步骤1．根据算法事先写出相应程序。

2．启动PC 机，进入tc 集成环境，输入代码。

3．编译调试。

4．调试通过，计算出正确结果后。

实验二解三对角线性方程组的追赶法⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=b b b x x x a a a a a a a a a n n nn n n n n b X A M M 2122122221112111....................................实验学时：2实验类型：验证实验要求：必修一、实验目的二、实验内容三、实验组织远行要求统一进行实验，一人一组四、实验条件PC机，tc2.0，Internet网五、实验步骤a)根据算法事先写出相应程序。

Runge-Kutta法求解方程

3. Runge-Kutta法求解方程3.1 题目用Runge-Kutta 4阶算法对初值问题y’=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1、0.2。

注：此方程的精确解为：y=e-20x3.2 方法介绍通过用一些点上函数值f(x ,y)的适当线性组合，来替代Euler法中的f(x k ,y k)，从而使方法的阶数更高，这就是Runge-Kutta法的基本思想。

在本次上机实习中采用的是标准4阶Runge-Kutta法，其计算形式如下：y k+1=y k+16[k1+2k2+2k3+k4]k1=hf(x k,y k)k2=hf(x k+h2,y k+k12)k3=hf(x k+h2,y k+k22)k4=hf(x k+h,y k+k3)3.3 计算结果及分析利用FORTRAN语言编写程序program R_K（附录4），计算结果如下。

表2 不同步长计算结果对比从上面结果中可以看出取步长h=0.1时候结果可以很好的收敛，与准确解之间的误差也越来越小，但是取h=0.2时候结果就发生了发散。

说明该方法对步长的选择比较敏感。

3.4 本题小结本题说明了采用Runge-Kutta法时候一定要选取合适的步长，否则可能造成结果很大的差异。

所以我们在采用Runge-Kutta法解题的时候一定要注意情况而选择不同的步长，甚至要试用多个步长。

4. 总结总的来说，通过完成本作业报告，我收获颇丰。

1）通过这次上机实习我知道了数值分析在计算中的重要性，我从实践的角度更为深入地掌握了拉格朗日插值法、分段低次插值法、SOR法、Runge-Kutta 法等数值方法。

刚开始对这些问题还不是很了解，结果在编程中让我渐渐的明白了这些方法的计算步骤，对这些方法有了更深入的了解。

2）让我对MATLAB和FORTRAN这两门语言有了更深入的理解。

在学习过程中本专业很多计算方法要利用数值分析。

所以我们勤加练习，在其中也发现了很多问题。

北航数值分析计算实习题目一幂法反幂法求矩阵特征值

《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。

1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。

2)使用反幂法求λs ，其中需要解线性方程组。

因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。

（2）与140k λλμλ-5011=+k 最接近的特征值λik 。

通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。

（3）2cond(A)和det A 。

1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。

2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。

由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。

2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。

数值分析计算实习第一题

直接用定义: ��(��)2 = ‖��‖2‖��−1‖2
求 A 的条件数很繁琐,需要先进行化简：
首先：
由于 A 是对称矩阵，
‖��‖2 = ��max(��)
说明：
1. 在所用的算法中，凡是要给出精度水平的ε，都取 ��=10−12。
2. 选择算法的时候应使矩阵 A 的所有零元素都不存储。
3. 打印以下内容：
(1)算法设计方案和思路。
(2)全部源程序。
(3)特征值��1,��501,��,��(��=1,2,⋯,39)以及��(��)2, det��的值(采用 e 型输出实型数，并至少显示 12 位有效数字)。
λi[16] -2.533970311130E+00 λi[38] 8.648666065193E+00
λi[17] -2.003230769563E+00 λi[39] 9.254200344575E+00
λi[18] -1.503557611227E+00 cond(A)2 1.925204273903E+03
λi[19] -9.935586060080E-01 det(A) 2.772786141752E+118

数值计算方法上机实习题答案.doc

1．设I n 1 x ndx ，0 5 x（ 1）由递推公式 I n 5I n 11，从 I 0的几个近似值出发，计算I 20；n解：易得： I 0 ln6-ln5=0.1823, 程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为： I 20= -3.0666e+010（ 2）粗糙估计 I 20，用 I n 1 1I n 1 1 ，计算 I 0；5 5n0.0079 1 x 20 1 x 200.0095因为dx I 20dx 6 5所以取 I 20 1(0.0079 0.0095) 0.0087 2程序为： I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endII 0= 0.0083（ 3）分析结果的可靠性及产生此现象的原因(重点分析原因 )。

首先分析两种递推式的误差；设第一递推式中开始时的误差为E0 I 0 I 0，递推过程的舍入误差不计。

并记 E n I n I n，则有 E n 5E n 1 ( 5) n E0。

因为 E20( 5) 20 E0 I 20，所此递推式不可靠。

而在第二种递推式中E0 1E1 (1)n E n，误差在缩小，5 5所以此递推式是可靠的。

出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。

2．求方程e x10x 2 0 的近似根，要求x k 1x k 5 10 4，并比较计算量。

（1）在 [0， 1]上用二分法；程序： a=0;b=1.0;while abs(b-a)>5*1e-4c=(b+a)/2;if exp(c)+10*c-2>0b=c;else a=c;endendc结果： c =0.0903（ 2）取初值x0 0，并用迭代 x k 1 2 e x ；10程序： x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果： x =0.0905（3）加速迭代的结果；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果： x =0.0995（ 4）取初值x00 ，并用牛顿迭代法；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x;end x 结果： x =0.0905（ 5）分析绝对误差。

数值分析计算实习题

1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-/*x^2+95549/72072*x-1/00*x^8-2168879/0*x^4+19/0*x^7+657859/*x^3+33983/ 0*x^5-13003/00*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x)输出结果为：S =/6464-2399/88*x+/1984*x^2+2656867/624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。

计算方法实验.

计算方法实验指导姓名______________学号______________院系______________专业______________哈尔滨工业大学计算方法实验指导根据实际问题建立的数学模型，一般不能求出所谓的解析解，必须针对数学模型的特点确定适当的计算方法，编制出计算机能够执行的计算程序，输入计算机，进行调试，完成运算，如果计算结果存在问题或不知是否正确，还需要重新确定新的计算方法，再编制出计算程序，输入计算机，重新调试，完成运算，直至获得正确的计算结果，这就是数值计算的全部过程。

学生在学习“计算方法”和“高级语言”等课程时普遍存在的问题是：只会套用教科书中的标准程序进行数值计算，很少有人能够独立地将学过的数值算法编制成计算机程序，至于灵活应用已经掌握的算法求解综合性较大的课题，则更是困难的事情。

编写《计算方法实验指导》的目的是：突出数值计算程序结构化的思想。

提高学生的编程能力，加深对“计算方法”课程内容的理解和掌握，为”计算方法“课程的教学服务，进一步奠定从事数值计算工作的基础。

具体地1.根据“计算方法”课程内容的特点，给出五个典型算法的分析流程，学生可以利用所掌握的“高级语言”顺利地编制出计算机程序，上机实习，完成实验环节的教学要求。

2.所有的计算实习题目都经过任课教师逐一检验，准确无误。

3.充分利用循环的思想、迭代的思想，给出算法结构描述和程序语言的对应关系，有利于学生编制相应的程序。

4.结合实习题目，提出实验要求，要求学生按规范格式写出相应的实验报告，实验报告成绩记入期末总成绩。

需要提醒学生：不能简单地套用现成的标准程序完成实验题目，应当把重点放在对算法的理解、程序的优化设计、上机调试和计算结果分析上，否则就失去实验课的目的啦。

5. 五个具体的实验题目是：实验题目1拉格朗日(Lagrange)插值实验题目2龙贝格(Romberg)积分法实验题目3四阶龙格—库塔(Runge—Kutta)方法实验题目4牛顿(Newton)迭代法实验题目5高斯(Gauss)列主元消去法要求必须完成其中三个（如果全部完成更好）。

计算方法实验报告

2019年计算方法（B）实验报告姓名：学号：专业：课程：计算方法（B）目录一、实验综述 (1)二、实验内容 (1)2.1 实验一 (1)2.2 实验二 (2)2.3 实验三 (3)2.4 实验四 (4)2.5 实验五 (6)三、思考总结 (7)附件A1 (8)附件A2 (9)附件A3 (10)附件A4 (12)附件A5 (14)一、实验综述计算方法在工程实践中得到了广泛的应用，是理工类研究生必备的知识技能。

按照2019年计算方法课程学习要求，本文对计算方法上机题目进行了算法设计、分析，利用matlab 2019b版本对算法进行实现，最终形成了实验报告。

以下为本次实验报告具体内容，包括五个实验部分和一个思考总结部分。

二、实验内容2.1 实验一2.1.1 实验题目用Jacobi迭代和Gauss-Seidel迭代解电流方程组，使各部分电流的误差均小于10-3。

2.1.2 算法分析a)首先列出方程组的系数矩阵A以及等式右端的矩阵b，A=[28,-3,0,0,0;-3,38,-10,0,-5;0,-10,25,-15,0;0,0,-15,45,0;0,-5,0,0,30 ];b=[10;0;0;0;0];为了验证A是否收敛，我们通过判断系数矩阵A是否为严格对角占优矩阵进行确定。

如果是，则可以进行Jacobi迭代和Gauss-Seidel迭代（利用matlab程序验证后，证明了矩阵A为严格对角占优矩阵）；如果不是，则需要采用其他方法进行判断迭代是否收敛。

b)对矩阵A分裂成三部分，，其中D为A的对角矩阵，E为A的下三角矩阵的相反数，F为A的上三角矩阵的相反数。

c) Jacobi迭代。

取x得初始向量为x=[0;0;0;0;0]，利用迭代公式进行循环计算，当的无穷范数小于10-3，即，停止循环。

d) Gauss-Seidel迭代。

取x得初始向量为x=[0;0;0;0;0]，利用迭代公式进行循环计算，当的无穷范数小于10-3，即，停止循环。

数值分析计算实习题列主元高斯消去法解线性方程组

数值分析计算实习题第5章解线性方程组的直接方法【选题列主元高斯消去法解线性方程组。

书上的计算实习题1、2、3都要求用列主元高斯消去法解线性方程组，所以考虑写一个普适的程序来实现。

对于线性方程组Ax二b,程序允许用户从文件读入矩阵数据或直接在屏幕输入数据。

文件输入格式要求:(1)第一行为一个整数n (2<=n<=100)，表示矩阵阶数。

(2)第2~n+l行为矩阵A各行列的值。

(3)第n+2~n+n+2行为矩阵b各行的值。

屏幕输入：按提示输入各个数据。

输出：A. b、det(A).列主元高斯消去计算过程、解向量X。

【算法说明】设有线性方程组Ax=b,其中设A为非奇异矩阵。

方程组的增广矩阵为«12«21[Nb] =第1步（k=l ）：首先在A的第一列中选取绝对值最大的元素®I,作为第一步的主元素:«|| H0然后交换（A, b）的第1行与第I行元素，再进行消元计算。

设列主元素消去法已经完成第1步到第k・l步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A（k）x=b（k）4? …4;)…唸)•忒••輕■[A.b]T[A ⑹,b")] = ••■咲■■■■■* *■〃伏)・• - %■第k步计算如下: 对于 k=l, 2, •…，0-1（1）按列选主元：即确定t使（2）如果tHk,则交换［A, b］第t行与第k行元素。

（3）消元计算54* J 叫=一鱼(=^ + 1,…,H)% 吗 <-«y + 〃如伽 (fJ = R + l,…/)b- <-勺+加汝仇， (i = /c + l,…,《)消元乘数mik 满足:n (%-D 内)X1 < ------ -- ---- 9(j = « 一 1,«一2■…J)tk M 1,(,=斤 +1, •••,«)fet e(4)回代求解【程序】/*【普适列主元消去法解线性方程组】对于线性方程组：Ax=b 输入：［选择屏幕直接输入］1.A的行阶数n(l <= 11<= 100)2.A的值3.b的值［选择读取文件1文件名(和主程序同级文件夹下)输出：I.A2・b3・ det(A)4解向疑X */#inciude <stdio.h>#include <stdlibJi>#include <niath.h> double A[1051(J05LA_B[I05][105Lb[105].x[105];double det A:int n.mark = 1;〃读入数据void input(){int ij:char ch[20],name[ 100];HLE *f;printf("\iv-An是否从文件读取数据(Y/N):”); scanf(”％st&ch);if(ch[0] = Y II ch[0] = y)( prinif(“请输入文件名(包括扩展需)：")； scanf("%s".name):f = fop cn(namc「T');fscanfCf/'%d'\&n);ford = 0:i < n;i 卄) for(j = 0;j < nJ ++) fscanf(f,'*%ir\&A[i]|j))：for(i = 0:i < n;i卄)fscanf(「％F・&b[i]);else{prin氓”请输入A的阶数：”)；scanf( '%d %d\&n)： prinifC请输入A的值：”)；for(i = 0:i V n:i ++)for(j = 0;j < n:j ++) scanf("%lf\&A[i]U]);phnifC请输入b的值：”)；for(i = 0:i < n;i 卄)scanf(''%lf\&b[i])：〃讣算行列式的值double det{double s[105](105] jni m){int z.jkdouble b[105][105Klotal = 0.r; /*b[Nl[N]用于存放，在矩阵s[Nl[N冲元素s[0]的余子式灯if(m>2){for(z = 0：z V m:z++){for(j = 0;j V m ・ 1 j ++)for(k = 0:k vn卜l;k ++)if(k >= z)bUKk] = sU+l](k+l];elsebLi][k] = s[j+l][k];if(z % 2==0)r=s[0)⑵ * dcKb.m - 1)：/*递归调用制elser= (-1) * s[0](z] * det(b.m -1);total = total + r;else if(m == 2)total = s[0][0] *s(l][l]-s[0](l] *s[l][01：else if(m == 1)total =s[0][0];return total;// 输出A^llb和dcl(A)void ouipui_l(){ int i j;primlTAW);for(i = 0;i < n:i ++){ for(j = 0:j < n:j ++)p rintf('-%15.4f\A[i]lj]);prinif(W);prinlf(5b = \rf);for{i = 0;i < n:i ++)prinif("%15.4l\iV\b[i]):printf("\ndet(A) = %・4f\n”・dclA);//主il•算函数void couni_x(){int ij,k:int max; double tmpjnik;〃构造增广矩阵for{i = 0;i<n:i++){for。

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

目录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 分解法求解线性方程组计算原理当n x R ∈是任意给定的非零向量，n v R ∈是任意给定的单位向量，则存在初等反射阵2T H I uu =-,使得Hx v σ=，其中σ为常数，当取单位向量2x vu x v σσ-=-P P 时，由u 确定的矩阵H 必定满足Hx v σ=，所以在计算过程中取u 的值为上述值。

计算方法上机作业

计算方法上机作业上机实习题目1.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。

在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。

已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0 1 2 3 4 5 6深度9.01 8.96 7.96 7.97 8.02 9.05 10.1 3分点7 8 9 10 11 12 13深度11.1812.2613.2813.3212.6111.2910.22分点14 15 16 17 18 19 20深度9.15 7.90 7.95 8.86 9.81 10.810.93(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；（1）算法思想分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。

在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。

计算光缆长度时，用如下公式：20()L f x ds=⎰200(f x =⎰191(k kk f x +==∑⎰= 本题采取三次样条插值的方法，因为三次样条插值方法是一种具有较好“光滑性”的分段插值方法。

根据提供的数据，只用x,y 值，不包含导数值，因此采用第三类三次插值多项式进行插值编程。

设计算法如下： 1.For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒2.For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13. 101h x x ⇒-4.For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2ba c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(114.3 i i d M ⇒+165.000;;c M d M d n n ⇒⇒⇒λnn n b a b ⇒⇒⇒2;;20μ6. 1111,γμ⇒⇒d b7.获取M 的矩阵元素个数，存入m8.For m k ,,3,2 =8.1kk k l a ⇒-1/μ8.2 k k k k c l b μ⇒⋅-1- 8.3 kk k k l d γγ⇒⋅-1-9.m m m M ⇒μγ/10. For 1,,2,1 --=m m k 10.1 k k k k k M M c ⇒⋅-+μγ/)(1 11. 获取x 的元素个数存入s 12. k ⇒113. For 1,,2,1-=s i 13.1 ifix x ≤~ then k i ⇒；breakelse k i ⇒+114.xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒---yh x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---（3）源程序clear; clc;x=0:1:20; %产生从0到20含21个等分点的数组 X=0:0.2:20;y=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93]; %等分点位置的深度数据n=length(x); %等分点的数目 N=length(X);%% 求三次样条插值函数s(x)M=y;for k=2:3; %计算二阶差商并存放在M中 for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1));a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i)for k=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);for m=1:N;k=1;for i=2:n-1if X(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);x2=X(m)-x(k); s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6))*x2)/H;end%% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光缆长度为 L=');disp(L);figureplot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置X/测量点','fontsize',16); %标注坐标轴含义ylabel('深度Y/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;（4）结果与分析铺设海底光缆的曲线图如下图所示：拟合结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图。

计算方法插值法(一)

浮点数的IEEE标准：舍入误差
当x在计算机中要求超过p个尾数位时，引起舍入误差
e* = x* - x （绝对误差）
计算机中的估计值
x*
e* < 2- p * 2E e Machine precision
相对误差 x* - x < 2- p = e
x
更多浮点数的内容见：
https:///wiki/Floating_point#IEEE_754:_floating_point_in_modern_computers
助教：赵渊明
上机实习作业
１.提交时间：作业布置下来两周内（如无特殊情况，晚交的作业做零分处理。有特殊情况的，需要提前得到授课老师许可，一事一议）（一般是周三）。２.提交内容：书面报告、源代码、源代码流程图及运行结果截图３.源代码要求：简洁、清楚、有较好的注释（助教能够运行程序并复制结果）。４.完成作业要求：鼓励同学之间讨论、合作完成作业，但最终程序、报告需要自己独立完成。如和别的同学交流过，请在上交作业中列出一起合作交流过的同学名字。如发现上交作业雷同，雷同作业做零分处理。 5. 编程语言：尽量选择Fortran或C语言，不建议使用Matlab。
| x1 x0 | 很小时
x1
15
由直线两点式可知，通过A，B的直线方程为
y
y0
y1 y0 x1 x0
(x
x0 )
L1( x)
称为线性插值(n=1的情况)
表示为如下对称形式： L1(x) y0l0 (x) y1l1(x)
其中，
l0(x)
x x1 x0 x1
l1 ( x)
x x0 x1 x0
精度提高的条件： ➢ 插值点与节点靠近 ➢ 内插精度一般比外推高 ➢ 插值点适当多