哈工大计算方法实验报告
计算方法_实验报告
一、实验目的1. 理解并掌握计算方法的基本概念和原理；2. 学会使用计算方法解决实际问题；3. 提高编程能力和算法设计能力。
二、实验内容本次实验主要涉及以下内容：1. 线性方程组的求解；2. 多项式插值；3. 牛顿法求函数零点；4. 矩阵的特征值和特征向量求解。
三、实验环境1. 操作系统：Windows 102. 编程语言：Python3.83. 科学计算库：NumPy、SciPy四、实验步骤及结果分析1. 线性方程组的求解（1）实验步骤a. 导入NumPy库；b. 定义系数矩阵A和增广矩阵b；c. 使用NumPy的linalg.solve()函数求解线性方程组。
（2）实验结果设系数矩阵A和增广矩阵b如下：A = [[2, 1], [1, 2]]b = [3, 2]解得：x = [1, 1]2. 多项式插值（1）实验步骤a. 导入NumPy库；b. 定义插值点x和对应的函数值y；c. 使用NumPy的polyfit()函数进行多项式拟合；d. 使用poly1d()函数创建多项式对象；e. 使用多项式对象计算插值点对应的函数值。
（2）实验结果设插值点x和对应的函数值y如下：x = [1, 2, 3, 4, 5]y = [1, 4, 9, 16, 25]拟合得到的二次多项式为：f(x) = x^2 + 1在x = 3时，插值得到的函数值为f(3) = 10。
3. 牛顿法求函数零点（1）实验步骤a. 导入NumPy库；b. 定义函数f(x)和导数f'(x)；c. 设置初始值x0；d. 使用牛顿迭代公式进行迭代计算；e. 判断迭代结果是否满足精度要求。
（2）实验结果设函数f(x) = x^2 - 2x - 3，初始值x0 = 1。
经过6次迭代，得到函数零点x ≈ 3。
4. 矩阵的特征值和特征向量求解（1）实验步骤a. 导入NumPy库；b. 定义系数矩阵A；c. 使用NumPy的linalg.eig()函数求解特征值和特征向量。
大学计算方法实验报告
《计算方法》实验报告实验题目实验报告1：非线性方程组的求解···················P1~2实验报告2：线性方程组解法·······················P3~4 实验报告3：Lagrange 插值多项式··················P5~7姓名：学号：班级：指导老师：时间：专业 序号 日期实验报告1：非线性方程组的求解【实验目的】1．用MATLAB 来实践进行牛顿法的变形，即对牛顿法进行了修正，使其应用更为方便，掌握用MATLAB 运用割线法求解非线性方程组。
2．运用MATLAB 进行隐函数作图。
【实验内容】[方法] 设a,b 为迭代初值,求两点(a,f(a)) 与 (b,f(b)) 的连线(割线)与 x 轴的交点记为 c ，再把迭代初值换成 b,c,重复计算.[要求] 把下面程序复制为新的 M-文件,去掉开头的 %再把 '?' 部分改写正确就是一个完整的程序,找前面一个例子试算【解】在牛顿迭代公式中用差商代替导数。
带入初值（a,f(a)）,（b,f(b)）,两点的连线与x 轴的交点作为c ，再把迭代初值换为b ，c ，重复计算。
【计算机求解】以y= x-exp(-x)为例初值a=0，b=1，误差不超过1.0*10^(-5)进行计算。
哈工大计算化学报告
三、分子动力学研究最新进展之蛋白质折叠
虽然不同的蛋白质序列不同， 但折叠背后的物理化学原理却是一 样的。 1998年，UCSF的Yong Duan 和Peter Kollman在Cray 3E超级 计算机上开展了第一个微秒级的蛋 白质折叠分子动力学模拟，有力地 推动了分子动力学的并行计算研究， 堪称分子动力学历史上的一次壮举!
2009年，Shan等人通过分子 动力学模拟Asp381不同质子化状 态下的构象变化，提出Asp381的 质子化是DFG翻转的关键，翻转的 功能是协助ATP的结合和ADP的释 放。
三、分子动力学研究最新进展之化 学反应动力学
大连化物所杨学明和张东辉两位研究员 在四原子反应动力学方面取得重要的进展 化学反应微分截面是反应碰撞中反应散 射到特定空间角度的反应物之间的有效碰撞 面积，是化学动力学实验能测量的最为精细 的物理量, 也是研究化学反应机理最为重要 的物理量之一。 发现在低碰撞能下激发态氟原子的反应 性居然比基态氟高出很多，表明波恩-奥本海 默近似的图像对于F2+D2反应是不适用的
三、分子动力学研究最新进展之化 学反应动力学
三、分子动力学研究最新进展之算 法进展
上海药物所蒋华良等利用国产 超级计算机开发了分子动力学模拟 的并行算法。
在并行计算中, 计算机将把密
度 矩阵的任务分配到并行计算机的各 个处理器上进行处理, 然后再把各 个计算结果合并处理得到最终计算 结果。 该研究成果在膜蛋白，药物设 计，细胞信号通路等方面研究具有 重大意义。
三、分子动力学研究最新进展之蛋白质动态 变化
2009年，，Dror等人开展了9个0.5~2.0 微 秒 的 分 子 动 力 学 模 拟 （ ～ 100， 000 原 子），通过分子动力学模拟推测，β2AR非 活性结构在非活性状态下（结合拮抗剂与 否），此离子键处于形成和消失的动态平衡 中。
哈工大计算机组成技术实验二
实验二顺序结构程序设计实验目的1 学会编制顺序结构的汇编语言程序2 进一步掌握完整汇编语言程序的结构、调试方法3 掌握运算指令对标志位的影响学会算术指令的格式与用法实验说明和注意事项1 应该熟练掌握由源文件到可执行文件过程2 掌握源程序构成的数据段和代码段之间的关系3 在DEBUG下运行的程序，程序中的结束语句最好不要用JMP $1、查找运行前数据区的内容E:\ MASM >DEBUG 2-3.EXE↙━U0 ↙分配数的据段用DEBUG调试、运行.EXE文件━D 1434:0 ↙分配的数据段运行前数据区内容源程序数据段内容用DEBUG调试、运行.EXE文件2、查找程序运行前各寄存器内容E:\ MASM >DEBUG 2-3.EXE↙━U0 ↙运行前各寄存器内容━R ↙用DEBUG调试、运行.EXE文件3、查找程序运行后数据区内容E:\ MASM >DEBUG 2-3.EXE↙━R ↙━G↙分配的数据段━D 1434 :0↙运行后数据区内容分配的数据段验收界面1、验收程序运行前和运行后数据区的内容 ━ D 1434 :0↙ ━ G ↙━ D 1434 :0↙验收的程序应是调试完成，没有错误。
今后在实验报告中要求填写“运行前、后数据区的内容”， 均验收此界面， 如果在CT2000 运行只验收程序运行后数据区内容的界面运行后运行前验收说明1、实验5表格填写，程序用T 命令运行运行结束后用D 命令观察Y 值━ T ↙ ━ T ↙执行ADD 前执行ADD 后执行ADD 前执行ADD 后━ T ↙（可继续执行）执行T 执行预习:分支及循环结构程序设计所有专业必做：实验2 、实验4 电类专业还要做实验1非电类专业还要做实验3选做：实验5 、实验6实验程序提前编写，在实验课上进行调试。
哈工大_数学实验报告
数学实验报告实验一Matlab的使用1.上机实验各种数据输入方法：程序语句：a=[1 2 3;4 5 6 ;7,8,9] 程序语句：linspace(1,10,5) 等等…………计算结果：a = 计算结果：ans =1 2 34 5 6 1.0000 3.2500 5.5000 7.7500 10.00007 8 92.(1) (a)方法：(b) 方法：程序语句：程序语句：a=[-3 5 0 8;1 -8 2 -1;0 -5 9 3;-7 0 -4 5]; a=[-3 5 0 8;1 -8 2 -1;0 -5 9 3;-7 0 -4 5];b=[0;2;-1;6]; b=[0;2;-1;6];inv(a)*b a\b计算结果：计算结果：ans = ans =-0.6386 -0.6386-0.4210 -0.4210-0.3529 -0.35290.0237 0.0237(2) 4个矩阵的生成语句：矩阵a 的生成语句：e=eye(3,3); a=[e r;o s]r=rand(3,2); 验证语句：o=zeros(2,3); a^2s=diag([1,2]);%此为一个任取的2X2 矩阵b=[e r+r*s; o s^2]计算结果相同：ans =1.0000 0 0 1.9003 1.45790 1.0000 0 0.4623 2.67390 0 1.0000 1.2137 2.28630 0 0 1.0000 00 0 0 0 4.00003.生成多项式的语句：poly ([2,-3,1+2i,1-2i,0,-6])计算结果：ans = 1 5 -9 -1 72 -180 0 计算x＝0.8,-x=-1.2 之值的指令与结果：指令：polyval([1,5,-9,-1,72,-180,0],0.8) 结果：ans= -100.2179指令：polyval([1,5,-9,-1,72,-180,0],-1.2) 结果：ans= 293.29004.求a的指令与结果：指令：a=compan([1,0,-6,3,-8])结果：a =0 6 -3 81 0 0 00 1 0 00 0 1 0求a的特征值的指令与结果：roots(p)的指令与结果为：指令：eig(a) 指令：roots([1,0,-6,3,-8])结果：结果：ans = ans =-2.8374 -2.83742.4692 2.46920.1841 + 1.0526i 0.1841 + 1.0526i0.1841 - 1.0526i 0.1841 - 1.0526i结论：利用友元阵函数a=company(p) 和eig(a) 可以与roots(p)有相同的作用，结果相同。
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。
例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。
在解决问题时，我们需要更有效地利用我们的知识和资源。
在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。
哈工大 计算机仿真技术实验报告 实验3 利用数值积分算法的仿真实验
实验3 利用数值积分算法的仿真实验(一、实验目的1) 熟悉MATLAB 的工作环境；2) 掌握MATLAB 的 .M 文件编写规则，并在命令窗口调试和运行程序； 3) 掌握利用欧拉法、梯形法、二阶显式Adams 法及四阶龙格库塔法构建系统仿真模型的方法,并对仿真结果进行分析。
二、实验内容系统电路如图2.1所示。
电路元件参数：直流电压源，电阻，电感，电容。
电路元件初始值：电感电流，电容电压。
系统输出量为电容电压。
连续系统输出响应的解析解为：))/sin (cos 1()(ωωωa t t e U t u at s c ⨯+⨯-=-（2-1）其中，LRa 2= ，221⎪⎭⎫⎝⎛-=L R LC ω 。
)(t u c 图2.1 RLC 串联电路三、实验要求1）利用欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta法构建系统仿真模型，并求出离散系统的输出量响应曲线；2）对比分析利用欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta 法构建系统仿真模型的仿真精度与模型运行的稳定性问题；3）分别编写欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta 法的.m 函数文件，并存入磁盘中。
.m 函数文件要求输入参数为系统状态方程的系数矩阵、仿真时间及仿真步长。
编写.m 命令文件，在该命令文件中调用已经编写完成的上述.m 函数文件，完成仿真实验；4）利用subplot 和plot 函数将输出结果画在同一个窗口中，每个子图加上对应的标题。
四、实验原理在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta 法等。
欧拉法、梯形法和二阶显式Adams 法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta 法是利用Taylor 级数匹配原理构造的仿真算法。
对于线性系统，其状态方程表达式为：()()()()()()t t t t t t ⎧=+⎨=+⎩xAx Bu y Cx Du 00)(x x =t (4-1) 式（4-1）中，[]Tn t x t x t x )()()(21 =x 是系统的n 维状态向量,[]Tm t u t u t u t )()()()(21 =u 是系统的m 维输入向量，[]Tr t y t y t y t )()()()(21 =y 是系统的r 维输出向量。
哈工大机械设计电算实验(matlab2010)
机械设计电算实验一：普通V带传动设计内容和任务1、普通V带传动设计内容给定原始数据：传递的功率P，小带轮转速n，传动比i及工作条件。
设计内容：带型号，基准长度Ld，根数Z，传动中心距a，带轮基准直径dd1、dd2，带轮轮缘宽度B，初拉力F0，和压轴力FQ。
2、CAD任务：（1）编制V带传动设计程序框图。
（2）编制V带传动设计原程序。
（3）按习题或作业中数据运行程序，要求对每一组数据各按三种V带型号计算，对每一种带型号选三种小带轮直径进行计算并输出所有结果。
二、变量标识符三、程序框图四、源程序与其说明程序说明：本程序用Matlab2010b软件编制，主要针对机械设计大作业上的题型设计。
使用时只要打开m文件，并点击运行，按照提示进行即可。
首先输入原始数据，然后根据自己的需要选择带型，中心距即可得到设计结果，无需再查找资料，方便高效,计算过程如有错误会进行提示，并返回到输入处进行改正。
而且该程序可以直接计算下一带轮直径或者计算下一带型，比较方便。
源程序如下(先复制到记事本，再新建一个m文件，粘贴)clear all;disp('欢迎使用本程序，请输入V带传动设计的原始数据');p=input('电动机工作功率(kw) P=');n=input('电动机满载转速(r/min) nm=');i=input('第一级传动比 i1=');a=input('请输入最短工作工作年限 a年b班 a=');b=input(' b=');disp('是否反复起动、正反转频繁或工作条件恶劣');ka1=input('是请输入1，否请输入0。
请输入：');disp('原动机类型:');disp('I类原动机包括普通笼型交流电动机，同步电动机，');disp(' 直流电动机（并激），n>=600r/min的内燃机')disp('II类原动机包括交流电动机(双笼型、滑环式、单相、大转差率)，');disp(' 直流电动机（复激、串激），单缸发动机，n<=600r/min的内燃机')d1=input('请选择原动机的类型，输入1或2。
计算方法上机实验报告
. 《计算方法》上机实验报告班级：XXXXXX小组成员：XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师：XXX二〇一八年五月二十五日前言通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。
以下为本次上机实验报告，按照实验内容共分为六部分。
实验一：一、实验名称及题目： Newton 迭代法例2.7(P38):应用Newton 迭代法求在附近的数值解，并使其满足.二、解题思路：设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标)(')(0001x f x f x x -=,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标)(')(1112x f x f x x -=称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把)(')(1n n n n x f x f x x -=+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。
三、Matlab 程序代码：function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0);y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1;while abs(x1-x0)>=tol x0=x1;y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; endx=double(x1) K四、运行结果：实验二：一、实验名称及题目：Jacobi 迭代法例3.7（P74）：试利用Jacobi 迭代公式求解方程组要求数值解为方程组的精确解. 二、解题思路：首先将方程组中的系数矩阵A 分解成三部分，即：U D L A ++=，D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。
哈工大数值分析报告上机实验报告
实验报告一题目： Gauss 列主元消去法摘要：求解线性方程组地方法很多,主要分为直接法和间接法.本实验运用直接法地Guass 消去法,并采用选主元地方法对方程组进行求解.前言：（目地和意义）1. 学习Gauss 消去法地原理.2. 了解列主元地意义.3. 确定什么时候系数阵要选主元数学原理：由于一般线性方程在使用Gauss 消去法求解时,从求解地过程中可以看到,若)1(-k kk a =0,则必须进行行交换,才能使消去过程进行下去.有地时候即使≠-)1(k kk a 0,但是其绝对值非常小,由于机器舍入误差地影响,消去过程也会出现不稳定得现象,导致结果不正确.因此有必要进行列主元技术,以最大可能地消除这种现象.这一技术要寻找行r ,使得)1()1(max ||->-=k ik ki k rk a a 并将第r 行和第k 行地元素进行交换,以使得当前地)1(-k kk a 地数值比0要大地多.这种列主元地消去法地主要步骤如下：1. 消元过程对k =1,2,…,n -1,进行如下步骤.1) 选主元,记ik ki rk a a >=max || 若||rk a 很小,这说明方程地系数矩阵严重病态,给出警告,提示结果可能不对.2) 交换增广阵A 地r ,k 两行地元素.kj rj a a ↔ (j=k,…,n +1)3) 计算消元kk kj ik ij ij a a a a a /-= (i=k+1,…,n ; j =k +1,……,n +1)2. 回代过程对k = n , n -1,…,1,进行如下计算)/(11,∑-=+-=nk j kk j kj n k k a x a a x至此,完成了整个方程组地求解.程序设计：本实验采用Matlab地M文件编写.Gauss消去法源程序：cleara=input('输入系数阵：>>\n')b=input('输入列阵b：>>\n')n=length(b);A=[a b]x=zeros(n,1);%%%函数主体for k=1:n-1;%%%是否进行主元选取if abs(A(k,k))<yipusilong;%事先给定地认为有必要选主元地小数yzhuyuan=1;else yzhuyuan=0;endif yzhuyuan;%%%%选主元t=A(k,k);for r=k+1:n;if abs(A(r,k))>abs(t)p=r;else p=k;endend%%%交换元素if p~=k;for q=k:n+1;s=A(k,q);A(k,q)=A(p,q);A(p,q)=s;endendend%%%判断系数矩阵是否奇异或病态非常严重if abs(A(k,k))< yipusilongdisp(‘矩阵奇异,解可能不正确’)end%%%%计算消元,得三角阵for r=k+1:n;m=A(r,k)/A(k,k);for q=k:n+1;A(r,q)=A(r,q)-A(k,q)*m;endendend%%%%求解xx(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1;s=0;for r=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s)x(k)=(A(k,n+1)-s)/A(k,k)end结果分析和讨论：例：求解方程⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡321643.5072.12-623.4712.31-32108-z y x .求解地结果为：x =[]367257386.0,05088607.0-49105822.0-， 例：求解方程⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡73109104-10172-42-4z y x 求得地结果为：x =[]857142857.1,89285714.0-196428571.0， 结论：采用Gauss 消去法时,如果在消元时对角线上地元素始终较大（假如大于10-5）,那么本方法不需要进行列主元计算,计算结果一般就可以达到要求,否则必须进行列主元这一步,以减少机器误差带来地影响,使方法得出地结果正确.实验报告二题目： Rung 现象产生和克服摘要：由于高次多项式插值不收敛,会产生Runge 现象,本实验在给出具体地实例后,采用分段线性插值和三次样条插值地方法有效地克服了这一现象,而且还取地很好地插值效果.前言：（目地和意义）1. 深刻认识多项式插值地缺点.2. 明确插值地不收敛性怎样克服.3. 明确精度与节点和插值方法地关系.数学原理：在给定n+1个节点和相应地函数值以后构造n 次地Lagrange 插值多项式,实验结果表明（见后面地图）这种多项式并不是随着次数地升高对函数地逼近越来越好,这种现象就是Rung 现象.解决Rung 现象地方法通常有分段线性插值、三次样条插值等方法.分段线性插值：设在区间[a, b ]上,给定n+1个插值节点a=x 0<x 1<…<x n =b和相应地函数值y 0,y 1,…,y n ,,求作一个插值函数)(x φ,具有如下性质：1) j j y x =)(φ,j=0,1,…,n .2) )(x φ在每个区间[x i , x j ]上是线性连续函数.则插值函数)(x φ称为区间[a, b ]上对应n 个数据点地分段线性插值函数.三次样条插值：给定区间[a, b ]一个分划⊿：a=x 0<x 1<…<x N =b若函数S(x)满足下列条件：1) S(x)在每个区间[x i , x j ]上是不高于3次地多项式.2) S(x)及其2阶导数在[a, b ]上连续.则称S(x)使关于分划⊿地三次样条函数. 程序设计流程：本实验采用Matlab 地M 文件编写.其中待插值地方程写成function 地方式,如下function y=f(x);y=1/（1+25*x*x ）;写成如上形式即可,下面给出主程序Lagrange 插值源程序：n=input('将区间分为地等份数输入：\n');s=[-1+2/n*[0:n]];%%%给定地定点,Rf为给定地函数x=-1:0.01:1;f=0;for q=1:n+1;l=1;%求插值基函数for k=1:n+1;if k~=q;l=l.*(x-s(k))./(s(q)-s(k));elsel=l;endendf=f+Rf(s(q))*l;%求插值函数endplot(x,f,'r')%作出插值函数曲线grid onhold on分段线性插值源程序clearn=input('将区间分为地等份数输入：\n');s=[-1+2/n*[0:n]];%%%给定地定点,Rf为给定地函数m=0;hh=0.001;for x=-1:hh:1;ff=0;for k=1:n+1;%%%求插值基函数switch kcase 1if x<=s(2);l=(x-s(2))./(s(1)-s(2));elsel=0;endcase n+1if x>s(n);l=(x-s(n))./(s(n+1)-s(n));elsel=0;endotherwiseif x>=s(k-1)&x<=s(k);l=(x-s(k-1))./(s(k)-s(k-1));else if x>=s(k)&x<=s(k+1);l=(x-s(k+1))./(s(k)-s(k+1));elsel=0;endendendff=ff+Rf(s(k))*l;%%求插值函数值endm=m+1;f(m)=ff;end%%%作出曲线x=-1:hh:1;plot(x,f,'r');grid onhold on三次样条插值源程序：（采用第一边界条件）clearn=input('将区间分为地等份数输入：\n');%%%插值区间a=-1;b=1;hh=0.001;%画图地步长s=[a+(b-a)/n*[0:n]];%%%给定地定点,Rf为给定地函数%%%%第一边界条件Rf"(-1),Rf"(1)v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4;for k=1:n;%取出节点间距h(k)=s(k+1)-s(k);endfor k=1:n-1;%求出系数向量lamuda,miula(k)=h(k+1)/(h(k+1)+h(k));miu(k)=1-la(k);end%%%%赋值系数矩阵Afor k=1:n-1;for p=1:n-1;switch pcase kA(k,p)=2;case k-1A(k,p)=miu(p+1);case k+1A(k,p)=la(p-1);otherwiseA(k,p)=0;endendend%%%%求出d阵for k=1:n-1;switch kcase 1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-miu(k)*v;case n-1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-la(k)*v;otherwised(k)=6*f2c([s(k) s(k+1) s(k+2)]);endend%%%%求解M阵M=A\d';M=[v;M;v];%%%%m=0;f=0;for x=a:hh:b;if x==a;p=1;elsep=ceil((x-s(1))/((b-a)/n));endff1=0;ff2=0;ff3=0;ff4=0;m=m+1;ff1=1/h(p)*(s(p+1)-x)^3*M(p)/6;ff2=1/h(p)*(x-s(p))^3*M(p+1)/6;ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p));ff4=Rf(s(p))-M(p)*h(p)*h(p)/6;f(m)=ff1+ff2+ff3+ff4 ;end%%%作出插值图形x=a:hh:b;plot(x,f,'k')hold ongrid on结果分析和讨论： 本实验采用函数22511)(xx f +=进行数值插值,插值区间为[-1,1],给定节点为 x j =-1+jh ,h=0.1,j =0,…,n .下面分别给出Lagrang e 插值,三次样条插值,线性插值地函数曲线和数据表.图中只标出Lagrang e 插值地十次多项式地曲线,其它曲线没有标出,从数据表中可以看出具体地误差.表中,L10(x)为Lagrang e插值地10次多项式,S10(x),S40(x)分别代表n=10,40地三次样条插值函数,X10(x),X40(x)分别代表n=10,40地线性分段插值函数.x f(x)L10(x)S10(x) S40(x) X10(x) X40(x) -1.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 -0.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.04244031830239 -0.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.04705882352941 -0.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.05245901639344 -0.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 -0.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.06639004149378 -0.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.07547169811321 -0.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.08648648648649 -0.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 -0.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.11678832116788 -0.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.13793103448276 -0.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.16494845360825 -0.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 -0.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.24615384615385 -0.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.35000000000000 0.30769230769231 -0.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.39024390243902 -0.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 -0.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.64000000000000 -0.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.80000000000000 -0.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.941176470588240 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.000000000000001.000000000000000.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.941176470588240.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.800000000000000.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.640000000000000.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.500000000000000.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.390243902439020.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.350000000000000.307692307692310.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.246153846153850.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.200000000000000.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.164948453608250.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.137931034482760.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.116788321167880.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.100000000000000.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.086486486486490.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.075471698113210.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.066390041493780.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.058823529411760.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.052459016393440.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.047058823529410.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.042440318302391.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154从以上结果可以看到,用三次样条插值和线性分段插值,不会出现多项式插值是出现地Runge现象,插值效果明显提高.进一步说,为了提高插值精度,用三次样条插值和线性分段插值是可以增加插值节点地办法来满足要求,而用多项式插值函数时,节点数地增加必然会使多项式地次数增加,这样会引起数值不稳定,所以说这两种插值要比多项式插值好地多.而且在给定节点数地条件下,三次样条插值地精度要优于线性分段插值,曲线地光滑性也要好一些.实验报告三题目： 多项式最小二乘法摘要：对于具体实验时,通常不是先给出函数地解析式,再进行实验,而是通过实验地观察和测量给出离散地一些点,再来求出具体地函数解析式.又因为测量误差地存在,实际真实地解析式曲线并不一定通过测量给出地所有点.最小二乘法是求解这一问题地很好地方法,本实验运用这一方法实现对给定数据地拟合. 前言：（目地和意义）1. 学习使用最小二成法地原理2. 了解法方程地特性 数学原理：对于给定地测量数据(x i ,f i )(i=1,2,…,n ),设函数分布为∑==mj j j x a x y 0)()(ϕ特别地,取)(x j ϕ为多项式j j x x =)(ϕ (j=0, 1,…,m )则根据最小二乘法原理,可以构造泛函∑∑==-=n i mj i j j i m x a f a a a H 110))((),,,(ϕ令0=∂∂ka H(k=0, 1,…,m ) 则可以得到法方程⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ求该解方程组,则可以得到解m a a a ,,,10 ,因此可得到数据地最小二乘解∑=≈mj j j x a x f 0)()(ϕ程序设计：本实验采用Matlab 地M 文件编写.其中多项式函数j j x =ϕ写成function 地方式,如下function y=fai(x,j)y=1;for i=1:jy=x.*y;end写成如上形式即可,下面给出主程序.多项式最小二乘法源程序clear%%%给定测量数据点(s,f)s=[3 4 5 6 7 8 9];f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05];%%%计算给定地数据点地数目n=length(f);%%%给定需要拟合地数据地最高次多项式地次数m=10;%%%程序主体for k=0:m;g=zeros(1,m+1);for j=0:m;t=0;for i=1:n;%计算内积(fai(si),fai(si))t=t+fai(s(i),j)*fai(s(i),k);endg(j+1)=t;endA(k+1,:)=g;%法方程地系数矩阵t=0;for i=1:n;%计算内积(f(si),fai(si))t=t+f(i)*fai(s(i),k);endb(k+1,1)=t;enda=A\b%求出多项式系数x=[s(1):0.01:s(n)]';y=0;for i=0:m;y=y+a(i+1)*fai(x,i);endplot(x,y)%作出拟合成地多项式地曲线grid onhold onplot(s,f,'rx') %在上图中标记给定地点结果分析和讨论：例用最小二乘法处理下面地实验数据.并作出)f地近似分布图.(x分别采用一次,二次和五次多项式来拟合数据得到相应地拟合多项式为：y1=-0.38643+0.82750x；y2=-1.03024+1.06893x-0.02012x2；y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5；分别作出它们地曲线图,图中点划线为y1曲线,实线为y2曲线,虚线为y5曲线.’x’为给定地数据点.从图中可以看出并不是多项式次数越高越好,次数高了,曲线越能给定点处和实际吻合,但别地地方就很差了.因此,本例选用一次和两次地多项式拟合应该就可以了.实验报告四题目： Romberg 积分法摘要：对于实际地工程积分问题,很难应用Newton-Leibnitz 公式去求解.因此应用数值方法进行求解积分问题已经有着很广泛地应用,本文基于Romberg 积分法来解决一类积分问题.前言：（目地和意义）1. 理解和掌握Romberg 积分法地原理；2. 学会使用Romberg 积分法；3. 明确Romberg 积分法地收敛速度及应用时容易出现地问题. 数学原理：考虑积分⎰=ba dx x f f I )()(,欲求其近似值,通常有复化地梯形公式、Simpsion公式和Cotes 公式.但是给定一个精度,这些公式达到要求地速度很缓慢.如何提高收敛速度,自然是人们极为关心地课题.为此,记T 1,k 为将区间[a,b ]进行2k 等分地复化地梯形公式计算结果,记T 2,k 为将区间[a,b ]进行2k 等分地复化地Simpsion 公式计算结果,记T 3,k 为将区间[a,b ]进行2k 等分地复化地Cotes 公式计算结果.根据Richardson 外推加速方法,可以得到收敛速度较快地Romberg 积分法.其具体地计算公式为： 1. 准备初值,计算)]()([21,1b f a f ba T +-=2. 按梯形公式地递推关系,计算∑-=-+-+-+-+=1201,11,11))5.0(2(221k i k k k k i ab a f a b T T 3. 按Romberg 积分公式计算加速值1441,11,11,--=----+---m mk m m k m m m k m T T T m=2,…,k4. 精度控制.对给定地精度R ,若R T T m m <--1,11,则终止计算,并取1,m T 为所求结果；否则返回2重复计算,直至满足要求地精度为止. 程序设计：本实验采用Matlab 地M 文件编写.其中待积分地函数写成function 地方式,例如如下function yy=f(x,y); yy=x.^3;写成如上形式即可,下面给出主程序Romberg 积分法源程序%%% Romberg 积分法 clear%%%积分区间 b=3; a=1;%%%精度要求 R=1e-5;%%%应用梯形公式准备初值 T(1,1)=(b-a)*(f(b)+f(a))/2; T(1,2)=T(1,1)/2+(b-a)/2*f((b+a)/2); T(2,1)=(4*T(1,2)-T(1,1))/(4-1); j=2; m=2;%%%主程序体%%%while(abs(T(m,1)-T(m-1,1))>R);%%%精度控制 j=j+1; s=0;for p=1:2^(j-2);s=s+f(a+(2*p-1)*h/(2^(j-1))); endT(1,j)=T(1,j-1)/2+h*s/(2^(j-1)); %%%梯形公式应用 for m=2:j; k=(j-m+1);T(m,k)=((4^(m-1))*T(m-1,k+1)-T(m-1,k))/(4^(m-1)-1); end end%%%给出 Romberg 积分法地函数表 I=T(m,1)结果分析和讨论： 1. 求积分dx x10063.精确解I= 24999676.运行程序得Romberg 积分法地函数表为1.0e+007 *4.70101520000000 3.05022950000000 2.63753307500000 2.49996760000000 2.49996760000000 0 2.49996760000000 0 0由函数表知Romberg 积分给出地结果为2.4999676*10^7,与精确没有误差,精度很高.2. 求积分dx xx⎰10sin . 直接按前面方法进行积分,会发现系统报错,出现了0为除数地现象.出现这种情况地原因就是当x=0时,被积函数分母出现0,如果用一个适当地小数ε（最好不要小于程序给定地最小误差值,但不能小于机器地最大精度）来代替可以避免这个问题.本实验取R =ε,可得函数表为：0.92073548319659 0.93979327500190 0.94451351171417 0.94569085359489 0.94598501993743 0.94614587227034 0.94608692395160 0.94608330088846 0.94608307538495 0 0.94608299406368 0.94608305935092 0.94608306035138 0 0 0.94608306038722 0.94608306036726 0 0 0 0.94608306036718 0 0 0 0故该函数地积分为0.94608306036718,取8位有效数字.3. 求积分dx x ⎰12sin本题地解析解很难给出,但运用Romberg 积分可以很容易给出近似解,函数表为：0.42073549240395 0.33406972582924 0.31597536075922 0.31168023948094 0.31062036680949 0.31035626065456 0.30518113697100 0.30994390573588 0.31024853238818 0.31026707591900 0.31026822526959 0 0.31026142365354 0.31026884083167 0.31026831215439 0.31026830189296 0 0 0.31026895856465 0.31026830376269 0.31026830173008 0 0 0 0.31026830119484 0.31026830172211 0 0 0 0 0.31026830172262 0 0 0 0 0故该函数地积分为0.31026830172262,取8位有效数字.结论：Romberg 积分通常要求被积函数在积分区间上没有奇点.如有奇点,且奇点为第一间断点,那么采用例3地方法,还是能够求出来地,否则,必须采用其它地积分方法.当然,Romberg 积分地收敛速度还是比较快地.。