计算方法上机作业
页眉内容 页脚内容1 计算方法上机报告
姓 名： 学 号： 班 级： 上课班级：页眉内容
页脚内容1 说明： 本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：
0681581482184161nnnnS
n，要求：
2. ① 若只需保留11个有效数字，该如何进行计算； 3. ② 若要保留30个有效数字，则又将如何进行计算；
（1） 算法思想 1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681nnnannnnn
；
2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数) （2）算法结构 1. ;0s
681581482184161nnnnt
n；
2. for 0,1,2,,ni if 10mt end; 3. for ,1,2,,0niii ;sst 页眉内容 页脚内容2 （3）Matlab源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0; for n=0:50 t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); if t<=10^(-m) %判断通项与精度的关系 break; end end; fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值 for i=n-1:-1:0 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s=s+t; %求和运算 end s=vpa(s,m) %控制s的精度
（4）结果与分析
当保留11位有效数字时，需要将n值加到n=7， s =3.1415926536； 当保留30位有效数字时，需要将n值加到n=22, s =3.14159265358979323846264338328。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 4. 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：
分点 0 1 2 3 4 5 6 深度 9.01 8.96 7.96 7.97 8.02 9.05 10.13 分点 7 8 9 10 11 12 13 页眉内容 页脚内容3 深度 11.18 12.26 13.28 13.32 12.61 11.29 10.22 分点 14 15 16 17 18 19 20 深度 9.15 7.90 7.95 8.86 9.81 10.80 10.93 ① 请用合适的曲线拟合所测数据点； ② 预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；
（1）算法思想 如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。 在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺
至B点，在某点处的深度为ih。
海底光缆线的长度预测模型 计算光缆长度时，用如下公式： 200()Lfxds
20'2
0()1()fxfxdx 页眉内容 页脚内容4 191'20()1()kkkfxfxdx
22()()xy
（2）算法结构 1. For ni,,2,1,0 1.1 iiMy
2. For 2,1k 2.1 For knni,,1, 2.1.1 ikiiiiMxxMM)/()(1
3. 101
hxx
4. For 1-,,2,1ni 4.1 11iii
hxx
4.2 bacchhhiiiiii2;1;)/(11
4.3 iidM1
6
5. 0000;;cMdMdnn
nnnbab2;;20 6. 1111
,db
7. 获取M的矩阵元素个数，存入m 8. For mk,,3,2 8.1 kkkla1
/
8.2 kkkkclb1
-
8.3 kkkkld1
-
9. mmm
M/ 页眉内容 页脚内容5 10. For 1,,2,1mmk 10.1 kkkkkMMc/)(1
11. 获取x的元素个数存入s 12. k1 13. For 1,,2,1si 13.1 if ixx~ then ki；break else ki1 14. xxxxxxhxxkkkkˆ~;~;11
yhxhMyxhMyxMxMkkkkkk~/]ˆ)6()6(6ˆ6[2211331
（3）Matlab源程序 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)); end end h(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); end M(1)=0; %选择自然边界条件 M(n)=0; 页眉内容 页脚内容6 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); end M(n)=y1(n)/u(n); %追赶法求解样条参数M(i) for k=n-1:-1:1; M(k)=(y1(k)-c(k)*M(k+1))/u(k); end s=zeros(1,N); for m=1:N; k=1; for i=2:n-1 if X(m)<=x(i); k=i-1; break; else k=i; end end H=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:N L=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2); end disp('所需光缆长度为 L='); disp(L); figure plot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图 xlabel('位置','fontsize',16); %标注坐标轴含义
线代上机实验报告(3篇)
第1篇一、实验目的1. 掌握线性代数基本概念和基本运算方法。
2. 熟悉MATLAB软件在解决线性代数问题中的应用。
3. 提高实际操作能力和编程能力。
二、实验环境1. 操作系统：Windows 102. 软件环境：MATLAB R2019b3. 实验设备：计算机三、实验内容1. 矩阵的基本运算2. 矩阵的秩3. 矩阵的逆4. 线性方程组的求解5. 特征值和特征向量6. 二次型及其标准形四、实验步骤1. 矩阵的基本运算（1）创建矩阵A：A = [1, 2, 3; 4, 5, 6; 7, 8, 9]（2）计算矩阵A的转置：A_transpose = A'（3）计算矩阵A的行列式：det_A = det(A)（4）计算矩阵A的逆：A_inverse = inv(A)2. 矩阵的秩（1）创建矩阵B：B = [1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12]（2）计算矩阵B的秩：rank_B = rank(B)3. 矩阵的逆（1）创建矩阵C：C = [1, 2; 3, 4]（2）判断矩阵C是否可逆：is_inverse = rank(C) == size(C, 1)（3）如果可逆，计算矩阵C的逆：C_inverse = inv(C)4. 线性方程组的求解（1）创建矩阵A和B：A = [1, 2; 3, 4]B = [5; 6]（2）使用MATLAB内置函数求解线性方程组：x = A \ B5. 特征值和特征向量（1）创建矩阵D：D = [4, 1; 2, 3]（2）计算矩阵D的特征值和特征向量：[V, D] = eig(D)6. 二次型及其标准形（1）创建矩阵E：E = [2, 1; 1, 3]（2）计算矩阵E的特征值和特征向量：[V, D] = eig(E)（3）将二次型E化为标准形：Q = V D inv(V)五、实验结果与分析1. 矩阵的基本运算（1）矩阵A：1 2 34 5 67 8 9（2）矩阵A的转置：1 4 72 5 83 6 9（3）矩阵A的行列式：（4）矩阵A的逆：-1.5 0.50.5 -0.52. 矩阵的秩矩阵B的秩为2。
上机报告-二分法,史蒂芬森迭代,割线法
计算方法上机实习报告[题目及目的要求]1.用二分法求方程0163=--x x 在[0,5]上根的近似值。
用牛顿迭代法求0133=--x x 在2=x 附近的实根。
2.完成史蒂芬森迭代加速法和割线法的子程序，并利用方程010423=--x x 对比对分法与一般迭代法。
[方法原理说明]1.二分法和牛顿迭代法：二分法是逐次把有根区间分半，舍弃无根区间而保留有根区间的一种逼近根的方法。
在这个过程中有根区间的长度以2的幂次方减少，当有根区间的长度小于给定的精度时，其中点就作为根的近似值。
牛顿迭代法的迭代格式为： 初值0x ()()k k k k x f x f x x '1-=+ （k=0,1,2....）显然，牛顿迭代格式能够迭代下去必须要求()k x f 的导数不能为0.当某个()0'=k x f 或很小时，迭代中断；当()k x f 满足一定条件时，牛顿迭代具有平方收敛速度。
该方法对初值0x 要求较高，若选取不当，则可能发散，若选取的好，则收敛很快。
2.史蒂芬森迭代加速法和割线法迭代法就是通过一个迭代格式进行反复迭代以产生一个序列。
若这个序列收敛于方程的根，就称这个迭代格式收敛。
史蒂芬森迭代加速法的迭代格式为()k k k k k k k x y z x y x x +---=+221()k k x f y = ，()k k y f z = （k=0,1,2....）割线法与牛顿迭代法一样，即在根的某个邻域内，()k x f 有直至二阶的连续导数，且()0'≠k x f ，则在邻域内选取初值10,x x ，迭代均收敛。
割线法的迭代格式为初值10,x x()()()()111--+---=k k k k k k k x x x f x f x f x x （k=2,3....）[计算步骤]1.二分法:1）给定a,b 及精度要求ep ； 2）计算x=（a+b ）/2 及()k x f ；3）若b-a<ep ，则返回主程序，x 作为近似根，否则转4； 4）若()()0<a f x f ，则b x ⇒，否则a x ⇒； 5）转2。
运筹学-大M法或两阶段法的上机实验
. 1实验报告实验课程名称运筹学实验工程名称大M法或两阶段法的上机实验年级专业学生学号00 学院实验时间：年月日实验容〔包括实验具体容、算法分析、源代码等等〕：1.书上P97页第6题：用大M 法和两阶段法求解以下线性规划问题。
ma* z=5；3213x x x ++ 约束条件：102x 4x x 321≥++，16.x 2x -x 321≤+A ：大M 法图1.1图1.2δ，得出目标函数的最优解*1=16，*2=0，由上面的结果可知，满足所求出的0≤j*3=0，s*4=16，R*5=0，s*=0，最优值是80。
当把M的值改为100000后，值还是一样的，这样就可以得出当M为100时，已经得出有效解。
B：两阶段法图1.3由图1.3可知，先进展线性规划的第一阶段，满足0≤j δ，且z 值为零，即说明存在一个可行解使得所有的人工变量都为零，此时*2=2.5，s*6=21，其余为0得出z=0。
接下来进展第二阶段，令z=5*1+*2+3*3-0s*4+0R*5+0s*6，和大M 的分析方法一样，最终将得到满足0≤j δ时到达最优解：当*1=16，*2=0，*3=0，s*4=6，R*5=0，s*6=0，最优值为80。
2.书上P97页第7题〔4〕大M 法和两阶段法求解以下线性规划问题 。
ma* z=；321x x 2x ++ 约束条件：，42x 2x 4x 321≥++，204x 2x 21≤+，162x 8x 4x 321≤++ A ：大M 法图2.1图2.2由上面的图 2.1可知，首先先输入数据即线性规划的系数如图 2.1所示令ma* z=321x x 2x ++-0s*4+0s*6+0s*7-MR*5；进展下一次迭代，以同样的方法一直下去，直到所求出的为止0≤j δ，就可以得出目标函数的最优解：*1=4，s*4=12，s*6=12，其余为0时，最优值为8。
当把M 的值改为100000后，值还是一样的，这样就可以得出当M 为100时，已经得出有效解。
数值分析2024上机实验报告
数值分析2024上机实验报告数值分析是计算数学的一个重要分支，它研究如何用数值方法来解决数学问题。
在数值分析的学习过程中，学生需要通过上机实验来巩固理论知识，并学会使用相应的数值方法来解决实际问题。
本篇报告将详细介绍2024年度数值分析上机实验的内容和结果。
一、实验内容2024年度数值分析上机实验分为四个部分，分别是：方程求根、插值与拟合、数值积分和常微分方程的数值解。
1.方程求根这部分实验要求使用数值方法求解给定的非线性方程的根。
常见的数值方法有二分法、牛顿法、割线法等。
在实验过程中，我们需要熟悉这些数值方法的原理和实现步骤，并对不同方法的收敛性进行分析和比较。
2.插值与拟合这部分实验要求使用插值和拟合方法对给定的一组数据进行拟合。
插值方法包括拉格朗日插值、牛顿插值等；拟合方法包括最小二乘拟合、多项式拟合等。
在实验中，我们需要熟悉插值和拟合方法的原理和实现步骤，并对不同方法的精度和稳定性进行比较。
3.数值积分这部分实验要求使用数值方法计算给定函数的积分。
常见的数值积分方法有梯形法则、辛普森法则、龙贝格积分等。
在实验过程中，我们需要熟悉这些数值积分方法的原理和实现步骤，并对不同方法的精度和效率进行比较。
4.常微分方程的数值解这部分实验要求使用数值方法求解给定的常微分方程初值问题。
常见的数值方法有欧拉法、改进的欧拉法、四阶龙格-库塔法等。
在实验中，我们需要熟悉这些数值解方法的原理和实现步骤，并对不同方法的精度和稳定性进行比较。
二、实验结果在完成2024年度数值分析上机实验后，我们得到了以下实验结果：1.方程求根我们实现了二分法、牛顿法和割线法，并对比了它们的收敛速度和稳定性。
结果表明，割线法的收敛速度最快，但在一些情况下可能会出现振荡；二分法和牛顿法的收敛速度相对较慢，但稳定性较好。
2.插值与拟合我们实现了拉格朗日插值和最小二乘拟合，并对比了它们的拟合效果和精度。
结果表明，拉格朗日插值在小区间上拟合效果较好，但在大区间上可能出现振荡；最小二乘拟合在整体上拟合效果较好，但可能出现过拟合。
数值分析大作业
数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。
上机作业4-Gaussian在结构化学中的应用二
昆明理工大学理学院上机实验报告课程名称：计算化学实验名称：专业班级：学生姓名：学号：上机作业4：Gaussian程序使用二：频率和热力学性质计算对乙烯酮( H2C=C=O)分子的振动频率和热力学性质进行计算。
1）写出乙烯酮( H2C=C=O)分子的内坐标及高斯输入文件。
，其中C=C: 1.35 C=O: 1.20 C-H: 1.09，H-C=C：120.0°。
（提示：设虚原子）%CHK=YIXITONG.CHK%rwf=yixitong.rwf#p B3LYP /6-31G* spYixitong0 1CC 1 R1X 2 R2 1 A1O 2 R3 3 A1 1 180.0H 1 R4 2 A2 3 0.0H 1 R4 2 A2 3 180.0R1=1.35R2=1.0R3=1.20R4=1.09A1=90.0A2=120.02）对H2C=C=O分子进行结构优化，给出结构的对称性，优化后的结构数据（键长、键角、二面角），能量值，并通过GaussView或ChemOffice将输入的结构图形以球棍形式列出。
计算方法：B3LYP基组设定：6-31G*优化：OPT对称性：C2V能量= -152.5984712921 C 0.0000002 C 1.314839 0.0000003 O 2.486220 1.171381 0.0000004 H 1.082743 2.077840 3.167280 0.0000005 H 1.082743 2.077840 3.167280 1.878582 0.0000003）对优化后的结构进行频率计算。
（提示：保持*.chk文件不变），找出计算的红外频率、振动模式及对称性，对频率进行矫正，并将校正后的数据和实验值相对比（填写表1）。
给出主要振动模式的图形、对应的峰值和红外光谱图（要求用origin作图）。
表1. 乙烯酮的红外分析振动模式实验值(cm-1) 计算值(cm-1) 校正值(cm-1) 红外强度对称性面内振动438 440.7177 423.662 2.9343 B2面外摇摆528 538.7049 517.8573 74.4313 B1面外摇摆588 590.7611 567.899 74.0799 B1面内扭曲977 1003.0527 964.22341 8.3269 B2伸缩振动1118 1179.8738 1134.2135 6.9638 A1面内伸缩1388 1434.7138 1379.1909 16.2242 A1面内拉伸2152 2239.9335 2153.2481 526.0268 A1面内拉伸3071 3209.1869 3084.9915 23.7501 A1面外拉伸3166 3298.0569 3170.4221 5.7145 B2计算方法：B3LYP基组设定：6-31G*频率计算：Freq geom=check guess=read频率矫正因子：0.9613(2)(1)(3)(4)(5) (6(7)(8)(9)-10010020030040050060070080005010015020025030010.7036549533175.0751177974257.3788746988E p s i l o nFrequency (cm-1)4）列表给出H2C=C=O分子的热力学数据。
计算方法与计算 实验一误差分析
% 输出的量--每次迭代次数k和迭代值xk,
%
--每次迭代的绝对误差juecha和相对误差xiangcha，
误差分析
误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算 中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。 因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时， 由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法 的好坏会影响到数值结果的精度。 一、实验目的
因为运行后输出结果为： y 1.370 762 168 154 49， yˆ =1.370 744 664 189
38， R 1.750 396 510 491 47e-005， WU= 1.782 679 830 970 664e-005 104 . 所
以， yˆ 的绝对误差为 10 4 ，故 y
③ 运行后输出计算结果列入表 1–1 和表 1-2 中。
④ 将算法 2 的 MATLAB 调用函数程序的函数分别用 y1=15-2*x^2 和
y1=x-(2*x^2+x-15)/(4*x+1)代替，得到算法 1 和算法 3 的调用函数程序，将其保
存，运行后将三种算法的前 8 个迭代值 x1, x2 ,, x8 列在一起（见表 1-1），进行
的精确解 x* 2.5 比较，观察误差的传播.
算法 1 将已知方程化为同解方程 x 15 2x2 .取初值 x0 2 ，按迭代公式
xk1 15 2xk2
C语言上机实验
实验一（第1章实验）实验目的：1.掌握运行C语言程序的全过程。
2.熟悉编译环境。
3.初步熟悉C语言程序的语法规定。
4.了解简单函数的使用方法。
实验内容：1.编程且上机运行：求3个整数的和。
2.编程且上机运行：求2个数的和、差、积和商。
3.编程且上机运行：输入3个数，求最大值。
4.编程且上机运行：输入圆的半径，求圆的面积和周长。
5.在屏幕上输出：“hello world!”实验结果：实验二（第3章实验）1.实验目的：理解C语言的类型系统。
实验内容：写程序测试数据－2在类型char，int，unsigned int，long int，unsigned long int 中存储情况。
实验过程：实验结果：参见各种类型的存储实现描述。
2.实验目的：了解混合类型计算中类型的转换规则。
实验内容：写程序测试多种类型数据一起运算时类型的转换及表达式结果的类型。
注意unsigned int和int数据运算时类型转换的方向。
实验过程：/** 类型转换问题* 试问下面两个表达式等价吗？*/#include <stdio.h>#include <stdlib.h>int main() {unsigned int ui,uj;ui = 1;uj = 2;if (ui < uj)printf("\n%u < %u is true !\n", ui, uj);elseprintf("\n%u < %u is false !\n", ui, uj);if (ui - uj < 0)printf("\n%u - %u <0 is true !\n", ui, uj);elseprintf("\n%u - %u <0 is false !\n", ui, uj);system("pause");return 0;}实验结果：参见类型转换规则。
计算方法上机程序
1.对分+扫描Private Function f(x!)f = x ^ 4 - 5 * x ^ 2 + x + 2 End FunctionPrivate Sub Form_Click() Dim a!, b!, h!, c!, p!, q!, x!a = InputBox("输入a")b = InputBox("输入b")h = InputBox("输入h")x = aDo While x < bIf f(x) * f(x + h) <= 0 Thenp = x: q = x + hDo While Abs(q - p) > 0.00001 c = (p + q) / 2If f(c) = 0 ThenExit DoElseIf f(p) * f(c) < 0 Thenq = cElsep = cEnd IfEnd IfLoopPrint "["; p, q; "]"; cEnd Ifx = x + hLoopEnd Sub2.用牛顿法求a的立方根，精度要求0.000005 Private Sub Form_Click()Dim x0 As Single, x1 As SingleDim a As Integera = InputBox("输入a")If a = 0 ThenPrint "a的立方根=0"EndEnd Ifx1 = (a ^ 1 / 3)Dox0 = (x1)x1 = (x0 - (x0 ^ 3 - a) / (3 * x0 ^ 2))Loop While (Abs(x1 - x0)) > 0.000005 Print "a的立方根为："; x1End Sub3.编写牛顿法求方程根1)x3-x2-2x-3=0(初值x0=2)Private Sub Form_Click()Dim x0 as Single,x1 as Singlex1=2Dox0=x1x1=x0-(x0^3-x0^2-2*x0-3)/(3*x0^2-2*x0-2) Loop while abs(x1-x0)>0.00001Print x1End sub2)x-sinx=1/2(初值x0=1)Private Sub Form_Click()Dim x0 As Single, x1 As Singlex1 = 1Dox0 = x1x1 = x0 - (x0 - Sin(x0) - 1 / 2) / (1 - Cos(x0))Loop While Abs(x1 - x0) > 0.00001Print x1End Sub4.列主元高斯消去法Private Sub Form_Click()Dim a(1 To 3, 1 To 4) As Single, t#, i!, j!, k!, r!, l#, x(1 To 3) As Single For i = 1 To 3For j = 1 To 4a(i, j) = InputBox("输入一个数")Print a(i, j);Next jPrintNext iFor k = 1 To 2r = kFor i = k + l To 3If Abs(a(i, k)) > Abs(a(r, k)) Then r = iNext iIf r <> k ThenFor i = 1 To 4t = a(k, i)a(k, i) = a(r, i)a(r, i) = tNext iEnd IfFor i = k + l To 3l = (a(i, k) / a(k, k))For j = k + l To 4a(i, j) = (a(i, j) - l * a(k, j)) Next jNext iNext kFor k = 3 To 1 Step -1s = 0For j = k + l To 3s = s + (a(k, j) * x(j)) Next jx(k) = (a(k, 4) - s) / a(k, k) Next kFor i = 1 To 3Print x(i),Next iEnd sub5.LU分解法Private Sub Form_Click()Const n = 4Dim a(1 To n, 1 To n) As Single, l(1 To n, 1 To n) As Single, u(1 To n, 1 To n) As SingleDim x(1 To n) As Single, y(1 To n) As Single, b(1 To n) As Single, s#, i!, j!, k!, r!For i = 1 To nFor j = 1 To na(i, j) = InputBox("输入a数组")Print a(i, j)Next jPrintNext iFor i = 1 To nb(i) = InputBox("输入b数组")Print b(i)PrintFor k = 1 To nFor j = k To ns = 0For r = 1 To k - 1s = s + l(k, r) * u(r, j)Next ru(k, j) = a(k, j) - sNext jFor i = k + 1 To ns = 0For r = 1 To k - 1s = s + (l(i, r) * u(r, k)) Next rl(i, k) = (a(i, k) - s) / u(k, k) Next iNext kFor i = 1 To ns = 0For k = 1 To i - 1s = s + l(i, k) * y(k)y(i) = b(i) - sNext iFor i = n To 1 Step -1s = 0For k = i + 1 To ns = s + (u(i, k) * x(k))Next kx(i) = (y(i) - s) / u(i, i)Next iFor i = 1 To nPrint x(i)Next iEnd Sub6.雅克比迭代Option Base 1Function cha(x!(), y!()) As SingleDim z As Single, i As Single, k As Singlen = 3z = Abs(x(1) - y(1))For i = 2 To nIf (z < Abs(y(i) - x(i))) Then z = Abs(x(i) - y(i))Next icha = zEnd FunctionPrivate Sub Form_Click()Dim a1, x(3) As Single, y(3) As SingleDim t As Single, s As Single, a(3, 3) As SingleDim i As Integer, j As Integer, k As Integer, n As Integer n = 3a1 = Array(10, -2, -1, -2, 10, -1, -1, -2, 5)b = Array(3, 15, 10)For i = 1 To n: y(i) = 0: Next ik = 1For i = 1 To 3For j = 1 To 3a(i, j) = a1(k)k = k + 1Next j, iFor k = 1 To 30For i = 1 To nx(i) = y(i)Next iFor i = 1 To nt = 0For j = 1 To nIf (i <> j )Then t = t + a(i, j) * x(j) Next jy(i) = ((b(i) - t) / a(i, i))Next iIf (cha(x, y) < 0.00001 )Then Print k;For i = 1 To nPrint y(i)Next iExit ForEnd IfNext kIf k > 30 Then Print "发散"End Sub7. 高斯－赛德尔迭代Option Base 1Function cha(x!(), y!()) As SingleDim z As Single, i As Single, k As Singlen = 3z = Abs(x(1) - y(1))For i = 2 To nIf z < Abs(y(i) - x(i)) Then z = Abs(x(i) - y(i))Next icha = zEnd FunctionPrivate Sub Form_Click()Dim a1, x(3) As Single, y(3) As SingleDim t As Single, s As Single, a(3, 3) As SingleDim i As Integer, j As Integer, k As Integer, n As Integer n = 3a1 = Array(10, -2, -1, -2, 10, -1, -1, -2, 5)b = Array(3, 15, 10)For i = 1 To n: x(i) = 0: Next ik = 1For i = 1 To 3For j = 1 To 3a(i, j) = a1(k)k = k + 1Next j, iFor k = 1 To 30For i = 1 To ny(i) = x(i)Next iFor i = 1 To nt = 0For j = 1 To nIf i <> j Then t = t + a(i, j) * x(j) Next jx(i) = (b(i) - t) / a(i, i)Next iIf cha(x, y) < 0.00001 Then Print k;For i = 1 To nPrint x(i)Next iEnd IfNext kIf k > 30 Then Print "发散"End Sub8.拉格朗日插值多项式，求在ｔ＝３.５处的函数值的近似值，节点由ｘ，ｙ数组给出Private Sub Form_Click()Const n = 3Dim p#, s!Dim x, y As Variantx = Array(1, 2, 3, 4)y = Array(4, 5, 14, 37)t = InputBox("input t ")p = 0For k = 0 To ns = 1For i = 0 To nIf (i <> k) Thens = s * ((t - x(i)) / (x(k) - x(i)))Next ip = p + (y(k) * s)Next kPrint pEnd Sub9.牛顿基本插值公式Private Sub Form_Click()Const n = 4Dim x(n) As Single, y(n) As Single, t#, p#, s# For i = 0 To (n)x(i) = InputBox("input x" & Trim(Str(i)))y(i) = InputBox("input y" & Trim(Str(i))) Next it = InputBox("input t")For k = 1 To nFor i = n To k(-１)y(i) = (y(i) - y(i - 1)) / ((x(i) - x(i - k)))Next ip = y(0)h = 1For i = 1 To nh = (h * (t - x(i - 1)))p = p + (h * y(i))Next iPrint "p="; pEnd Sub10.拟合Private Sub Form_Click()Dim l#, m#, n#, i%, j%, k%, t1#Dim x As Variant, y As Variantn = 7m = 2x = Array(0, 1, 2, 3, 4, 5, 6, 7)y = Array(0, 5, 3, 2, 1, 2, 4, 7)ReDim a(0 To m, 0 To m + 1) As Single, t(n) As Single For i = 0 To mFor k = 1 To ns = s + x(k) ^ i * y(k) Next ka(i, m + 1) = sFor j = 0 To ms = 0For k = 1 To ns = s + x(k) ^ (i + j) Next ka(i, j) = sNext jNext iFor i = 0 To mFor j = o To m + 1 Print a(i, j),Next jPrintNext iFor k = 0 To mr = kFor i = k + 1 To mIf Abs(a(i, k)) > Abs(a(r, k)) Then r = iEnd IfNext iIf r <> k ThenFor i = 0 To m + 1t1 = a(k, i)a(k, i) = a(r, i)a(r, i) = t1Next iEnd Ifl = 1For i = k + 1 To ml = a(i, k) / a(k, k)For j = k + 1 To m + 1a(i, j) = a(i, j) - l * a(k, j)Next jNext iNext kFor k = m To step - 1s = 0For j = k + 1 To ms = s + a(k, j) * t(j)Next jt(k) = (a(k, m + 1) - s) / a(k, k) Next kPrint "y="; t(0)For i = 1 To mIf t(i) >= 0 Then Print "+" Print t(i); "*x^;i;"Next iEnd Sub11.SimpsonFunction f!(x!)f = x + x * Exp(x)End FunctionPrivate Sub Form_Click() Dim b!, N!, h!, x!, s!, a!a = 0:b = 1: N = 8h = (b - a) / (2 * N)x = as = f(a)For i = 1 To Nx = x + hs = s + 4 * f(x)x = x + hs = s + 2 * f(x)Next is = h / 3 * (s - f(b))Print sEnd Sub（2）Private Sub form_click()Dim a As Single, b As Single, eps As Single, s As Single Dim x As Single, h As Single, t1 As Single, t As Singlea = InputBox("输入积分下限a")b = InputBox("输入积分上限b")eps = InputBox("输入精度要求eps")h = b - at2 = (h / 2) * (f(a) + f(b))Dot1 = t2s = 0For x = a + h / 2 To b Step hs = s + f(x)Next xt2 = t1 / 2 + (h / 2) * sh = h / 2Loop While Abs(t1 - t2) > eps Print "积分的近似值："; t2End SubFunction f(x As Single) As Singlef = 1 / (1 + x * x)End Function12．用牛顿法求方程的附近根Private Sub Form_Click()Dim x#, x1#x1 = 1Dox = x1x1 = x - (x - Exp(-x)) / (1 + Exp(-x)) Loop While Abs(x1 - x) > 10 ^ (-5) Print x1End Sub。
计算方法编程作业1_拉格朗日插值与牛顿插值
西华数学与计算机学院上机实践报告课程名称：计算方法年级：2012级上机实践成绩：指导教师：严常龙姓名：贺容英上机实践名称：拉格朗日插值和牛顿插值法学号：上机实践日期：yyyy.mm.dd 上机实践编号：1312012070102209 上机实践时间：2014.10.27一、目的1．通过本实验加深对拉格朗日插值和牛顿插值法构造过程的理解；2．能对上述两种插值法提出正确的算法描述编程实现。
二、内容与设计思想自选插值问题，编制一个程序，分别用拉格朗日插值法和牛顿插值法求解某点的函数近似值。
（从课件或教材习题中选题）已知y=f(三、使用环境操作系统：win7软件环境：vs2012四、核心代码及调试过程4.1核心代码1、拉格朗日插值法代码如下double lagrangesf(point points[],int t){int n=t;int i,j;double x,tmp=1,lagrange=0;printf("请输入需要计算的x的值：");scanf("%lf",&x);for(i=0;i<=n-1;i++){tmp=1;for(j=0;j<=n-1;j++){if(j!=i)tmp=tmp*(x-points[j].x)/(points[i].x-points[j].x);}lagrange=lagrange+tmp*points[i].y;}printf("lagrange(%lf)=%lf\n",x,lagrange);return 0;}2、牛顿插值法代码如下double newtonsf(point points[],int t){int n=t;int i,j;double d[maxt+1];double x,tmp,newton=0;printf("差商表\n");printf("***************************************************\n"); printf("x ");for(i=0;i<=n-1;i++){printf("%lf ",points[i].x);}printf("\n");printf("y ");for(i=0;i<=n-1;i++){d[i]=points[i].y;printf("%lf ",points [i].y);}printf("\n");for(i=0;i<n-1;i++){printf("%d阶差商",i+1);for(int t=1;t<=i+12;t++)printf(" ");for(j=n-1;j>i;j--){d[j]=(d[j]-d[j-1])/(points[j].x-points[j-i-1].x);//计算差商printf("%lf ",d[j]);}printf("\n");}printf("***************************************************\n"); printf("请输入需要计算的x的值：");scanf("%lf",&x);tmp=1;newton=d[0];for(i=0;i<n-1;i++){tmp=tmp*(x-points[i].x);newton=newton+tmp*d[i+1];}printf("newton(%lf)=%lf\n",x,newton);return 0;}3、主函数中负责输入被插值点的输入，以及拉格朗日插值法和牛顿插值法的调用，代码如下int n=0;int i,j;point points[maxt+1];double x,tmp=0,lagrange=0;do{printf("请输入被插值点数目：");scanf("%d",&n);if(n>maxt){printf("被插值点数超出范围%d",maxt);return 0;}}while(n<=0);printf("请输入被插值点:\n");for(i=0;i<=n-1;i++){scanf("%lf%lf",&points[i].x,&points[i].y);}printf("lagrange插值\n");lagrangesf(points,n);printf("newton插值\n");newtonsf(points,n);system("pause");4.2、调试过程1、在拉格朗日插值法调试过程中，累乘过程中用来承载累乘的tmp没有重新赋值为1，导致结果始终不正确。