数值计算方法上机实习题
数值分析计算实习题二
《数值分析》计算实习题二算法设计方案1．主要计算步骤：计算函数f(x,y)在拟合所需的节点处的函数值。
将各拟合节点（x i，y j）分别带入非线性方程组0.5 cos t + u + v + w – x = 2.67t + 0.5 sin u + v + w – y = 1.070.5t + u + cos v + w – x =3.74t + 0.5u + v + sin w – y =0.79解非线性方程组得解向量（t ij，u ij，v ij，w ij）。
对数表z(t,u)进行分片二次代数插值，求得对应（t ij，u ij）处的值，即为f(x i，y j) 的值。
对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。
每次拟合后验证误差大小，直到满足要求。
2．求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。
3．对z(t，u)进行插值选择分片二次插值。
4．拟合基函数φr(x)ψs(y)选择为φr(x)=x r，ψs(y)=y s。
拟合系数矩阵c通过连续两次解线性方程组求得。
一．源程序#include "stdio.h"#include "stdlib.h"#include "math.h"void Doolittle(double *A,int n,int *M)//功能说明：对n阶矩阵A进行选主元的Doolittle分解//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果// 存储于A中// n：方阵A的维数// M；（返回参数）n维向量，记录选主元过程中行交换的次序{int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double));for(k=0;k<n;k++){for(i=k;i<n;i++){s[i]=A[i*n+k];for(t=0;t<k;t++) s[i]-= A[i*n+t] * A[t*n+k];}Maxs=abs(s[k]); M[k]=k;for(i=k+1;i<n;i++){if(Maxs<abs(s[i])){Maxs=abs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;t<n;t++){temp=A[k*n+t];A[k*n+t]=A[M[k]*n+t];A[M[k]*n+t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;}if(Maxs<(1e-14)){s[k]=1e-14;printf("%.16e方阵奇异\n",Maxs);}A[k*n+k]=s[k];for(j=k+1;(j<n)&&(k<n-1);j++){for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j];A[j*n+k]=s[j]/A[k*n+k];}}}void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中//参数说明：A：经Doolittle分解后的方阵// n：方阵A的维数// b：方程组的右端向量// x：（返回参数）方程组的解向量{int i,t;for(i=0;i<n;i++){x[i]=b[i];for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t];}for(i=n-1;i>-1;i--){for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t];x[i]/=A[i*n+i];}}void Transpose(double *A,int m,int n,double* AT)//功能说明：求m×n阶矩阵A的转置AT//参数说明：A：已知m×n阶矩阵// m：A的行数// n：A的列数// AT：（返回参数）A的转置矩阵（n×m）{int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];}void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。
数值分析上机实习题
2019-2020 第1学期数值分析上机实习题总目标：会算，要有优化意识。
（以下程序要求以附件1例题代码格式给出）1. 对给定的线性方程组Ax b =进行迭代求解。
（1）给出Jacobi 迭代的通用程序。
（2）给出Gauss-Seidel 迭代的通用程序。
调用条件：系数矩阵A ，右端项b ，初值0x ，精度要求ε。
输出结果：方程组的近似解。
给定线性方程组211122241125x --⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪---⎝⎭⎝⎭，和122711122215x -⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，取初值0x 为0， 分别利用Jacobi 迭代和G-S 迭代进行求解，观察并解释其中的数学现象。
2. 利用紧凑格式（即直接分解法或逐框运算法）对给定的矩阵A 进行Doolittle 分解，并用其求线性方程组的解。
调用条件：矩阵A 。
输出结果：单位下三角矩阵L 和上三角矩阵U 。
给定矩阵1112A ⎛⎫= ⎪⎝⎭，利用以下算法：1）将A 作Doolittle 分解11A LU =，2）令211A U L =，并对2A 作Doolittle 分解222A L U =，3）重复2）的过程令11n n n A U L --=，并对n A 作Doolittle 分解n n n A L U =，2,3,4,n =， 观察n L ，n U ，n A 的变化趋势，思考其中的数学现象。
3. 给定函数21(),12511f x x x -≤+≤=，取164,8,n =，用等距节点21,i i n x =-+ 0,1,,1i n =+对原函数进行多项式插值和五次多项式拟合，试画出插值和拟合曲线，并给出数学解释。
4. 给出迭代法求非线性方程()0f x =的根的程序。
调用条件：迭代函数()x ϕ，初值0x输出结果：根的近似值k x 和迭代次数k给定方程32()10f x x x =--=，用迭代格式1k x +=0 1.5x =附近的根，要使计算结果具有四位有效数字，利用估计式*1||1||k k k L x x x x L -≤---，或估计式*10||1||kk L x x x x L-≤--来判断需要的迭代次数，分别需要迭代多少次？两者是否有冲突？5. 利用数值求积算法计算()ba f x dx ⎰。
数值计算方法上机
数值计算方法上机实习题1． 设⎰+=105dx xx I nn ， （1） 由递推公式n I I n n 151+-=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用nI I n n 51511+-=-，计算0I ；（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
解：（1）由题可得：0I =ln6-ln5=0.1823， 源程序为：>> I=0.182;>> for n=1:20, I=(-5)*I+1/n; end >> I运行结果为：I = -3.0666e+010 即20I =-3.0666e+010（2）由12000.00796x dx ≈⎰，1200.00955x dx ≈⎰1120202065x x dx I dx ⎰⎰20I =1(0.00790.0095)2+=0.0087源程序为：>> I=0.0087; >> for n=1:20,I=(-1/5)*I+1/(5*n); end >> I运行结果为：I = 0.0083 即0I =0.0083(3) 首先分析两种递推式的误差；设第一递推式中开始时的误差为'000E I I =-，递推过的舍入误差不计。
并记'n n n E I I =-，则有105(5)n n n E E E -=-==-。
因为2020020(5)E E I =-，所此递推式不可靠。
而在第二种递推式中0111()55n n E E E =-==-，误差在缩小，所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2． 求方程0210=-+x e x的近似根，要求41105-+⨯<-k k x x ，并比较计算量。
（1） 在[0，1]上用二分法； （2） 取初值00=x ，并用迭代1021x k e x -=+；（3） 加速迭代的结果；（4） 取初值00=x ，并用牛顿迭代法； （5） 分析绝对误差。
数值计算方法上机实习题考证
数值计算⽅法上机实习题考证--------------------------------------------------- 此⽂档包含我们计算⽅法的经典算法包含（数值计算⽅法上机实习题）1．设?+=105dx xx I nn ，（1）由递推公式n I I n n 151+-=-，从0I 的⼏个近似值出发，计算20I ；（2）粗糙估计20I ，⽤nI I n n 51511+-=-，计算0I ；（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。
(1) 解答：n=0，0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=x d xdx x dx x x I nn这⾥可以⽤for 循环，while 循环，根据个⼈喜好与习惯：for 循环程序： While 循环程序： I=0.1823; I=0.1823; for n=1:20 i=1;I=(-5)*I+1/n; while i<21 End I=(-5)*I+1/i; I i=i+1; fprintf('I20=%f',I) end I = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009 (2) 粗略估计I 20： Mathcad 计算结果： for 循环程序： While 循环程序： >> I=0.007998; I=0.007998; >> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> II =0.0083(3) 算法误差分析：计算在递推过程中传递截断误差和舍⼊误差第⼀种算法：（从1——>20）1x x 205x +d 7.998103-?=*000e I I =-*115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===误差放⼤了5n倍，算法稳定性很不好；第⼆种算法：（从20——>1）*n n ne I I =-***111111111()()555555n n n n n n nn e I I I I I I e n n ---=-=-+--+=-=0111...()55n ne e e ===误差在逐步缩⼩，算法趋近稳定，收敛。
数值计算上机作业
西南交通大学数值计算课程上机作业组号第1组组员学号班级任课教师2017年12月10日目录第1题........................................................................................................ 错误!未定义书签。
1.1 (a) ..................................................................................................... 错误!未定义书签。
1.2 (b) ..................................................................................................... 错误!未定义书签。
第2题. (1)2.1 (a)................................................................................................. 错误!未定义书签。
2.2 (b) ................................................................................................ 错误!未定义书签。
第3题. (1)第4题 (2)4.1 (a)................................................................................................. 错误!未定义书签。
4.2 (b) ................................................................................................ 错误!未定义书签。
数值计算上机实习题目（matlab编程）
数值计算上机实习题目（matlab编程）非线性方程求根一、实验目的本次实验通过上机实习，了解迭代法求解非线性方程数值解的过程和步骤。
二、实验要求1、用迭代法求方程230x x e -=的根。
要求：确定迭代函数?(x)，使得x=?(x)，并求一根。
提示：构造迭代函数2ln(3)x ?=。
2、对上面的方程用牛顿迭代计算。
3、用割线法求方程3()310f x x x =--=在02x =附近的根。
误差限为410-，取012, 1.9x x ==。
三、实验内容1、（1）首先编写迭代函数，记为iterate.mfunction y=iterate(x)x1=g(x); % x 为初始值。
n=1;while(abs(x1-x)>=1.0e-6)&(n<=1000) % 迭代终止的原则。
x=x1;x1=g(x);n=n+1;endx1 %近似根n %迭代步数（2）后编制函数文件?(x)，记为g.mfunction y=g(x)y=log(3*x.^2);（3）设初始值为0、3、-3、1000，观察初始值对求解的影响。
将结果记录在文档中。
>>iterate(0)>>iterate(3) 等等2、（1）首先编制牛顿迭代函数如下,记为newton.mfunction y=newton(x0)x1=x0-fc(x0)/df(x0); % 牛顿迭代格式n=1;while(abs(x1-x0)>=1.0e-6)&(n<=1000000) % 迭代终止的原则。
x0=x1;x1=x0-fc(x0)/df(x0);n=n+1;endx1 %近似根n %迭代步数(2)对题目中的方程编制函数文件，记为fc.mfunction y=fc(x)y=3*x.^2-exp(x)编制函数的导数文件，记为df.mfunction y=df(x)y=6*x-exp(x)(3)在MATLAB 命令窗计算，当设初始值为0时，newton(0)；给定不同的初始值，观察用牛顿法求解时所需要的迭代步数，并与上面第一题的迭代步数比较。
数值计算方法上机实验考试答案
数值计算方法上机实验考试答案1. 小型火箭初始质量为900千克，其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力。
当燃料用尽时引擎关闭。
设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）。
重力加速度取9.8米/秒2.A. 建立火箭升空过程的数学模型（微分方程）；B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度。
解：A ． 建立模型：)(t x ——t 时刻的火箭高度；T =30000（牛顿）——火箭推力，当时间t >40秒时，T=0； m t m 150-=——火箭飞行过程中的质量，t >40秒时，300=m 千克0m =900（千克）——火箭初始质量； 1m =600（千克）——燃料质量；c =15（千克/秒）——燃料的燃烧速率； k =0.4（千克/米）——空气阻力系数； g =9.8（米/秒2）——重力加速度由能量守恒定律，可得到火箭飞行过程的方程：mg T t x k t x m -+'-=''2)]([)(这是一个初值问题，初始条件为 0)0(,0)0(='=x x设)(),(21t x x t x x '==，则问题化为求下列微分方程组的初值问题：⎪⎩⎪⎨⎧--+--='='g t m T x t m k x x x 151******** 0)0(,0)0(21==x xB ． 关闭引擎时4015/600==t （秒），所求的是此时火箭的高度)40(1x ；速度)40(2x ；加速度)40()40(21x x '''或，及火箭到最高点的时间m t 和高度)(1m t x 。
具体的Matlab 程序如下： 首先建立微分方程的的m 文件：function y=huojian(t,x)k=0.4;g=9.8;m0=900;T=30000;m=m0-15*t;if t>40T=0;m=300;endy=[x(2),-(k/m)*x(2)^2+T/m-g]';主程序：%feixing.mk=0.4;g=9.8;m0=900;T=30000;x0=[0,0];ts=0:1:55;[t,x]=ode23('huojian',ts,x0);[t,x(:,1)]%------a=[t,x];x40=a(41,2) %燃料用尽时的高度v40=a(41,3) %燃料用尽时的速度a40=-(k/300)*v40^2+T/300-g %燃料用尽时的加速度%-------xmax=max(x(:,1)) %火箭到达最高点的高度subplot(2,1,1),plot(t,x(:,1)),title('altitude') subplot(2,1,2),plot(t,x(:,2)),title('speed')运行结果为：1.0e+003 *0 00.0010 0.01180.0020 0.04750.0030 0.10670.0040 0.18890.0050 0.29270.0060 0.41680.0070 0.55920.0080 0.71800.0090 0.89140.0100 1.07700.0380 7.80510.0390 8.06310.0400 8.32240.0410 8.54000.0420 8.70040.0430 8.82420.0440 8.92180.0450 8.99940.0460 9.06070.0470 9.10830.0480 9.14370.0490 9.16810.0500 9.18220.0510 9.18640.0520 9.18070.0530 9.16510.0540 9.13900.0550 9.1018x40 =8322.4v40 = 254.1728a40 = 4.0616xmax =9186.4关闭引擎时4015/600==t （秒），此时火箭的高度h=)40(1x =8322.4米，速度v=)40(2x =254.1728米/秒，加速度为a=)40(1x ''= 4.0616米/秒2，火箭到最高点的时间m t =51秒，高度)(1m t x =9186.4米。
计算方法上机实习题大作业(实验报告)
计算方法实验报告班级： 学号： 姓名： 成绩：1 舍入误差及稳定性一、实验目的（1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；（2）通过上机计算，了解舍入误差所引起的数值不稳定性二、实验内容1、用两种不同的顺序计算1000021n n -=∑，分析其误差的变化 2、已知连分数()101223//(.../)n n a f b b a b a a b =++++，利用下面的算法计算f ： 11,i n n i i i a d b d b d ++==+ (1,2,...,0)i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f3、给出一个有效的算法和一个无效的算法计算积分1041nn x y dx x =+⎰ (0,1,...,10)n = 4、设2211N N j S j ==-∑，已知其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭（1）编制按从大到小的顺序计算N S 的程序（2）编制按从小到大的顺序计算N S 的程序（3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数三、实验步骤、程序设计、实验结果及分析1、用两种不同的顺序计算1000021n n -=∑，分析其误差的变化 （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.逆序计算结果分析：两种不同顺序计算结果是一样的，顺序计算误差从一开始就很小，而逆序计算误差最开始十分大，后来结果正确。
计算方法上机实习题
数值计算方法上机实习题1． 设⎰+=105dx xx I nn ， （1） 由递推公式nI I n n 151+-=-，从0=0.1822I , 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用nI I n n 51511+-=-，计算0I ；（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
解：（1）程序如下: clear all clc I=0.1822; %题中的已知数据 for n=1:20; I=(-5)*I+1/n; %由递推公式所得 end fprintf('I20=%f\n',I) M=0.1823; %与I 的计算结果形成对比for i=1:20; M=(-5)*M+1/i; %由递推公式所得 end fprintf('M20=%f\n',M) 输出结果为： I20=-11592559237.912731 M20=-2055816073.851284 （2）程序如下: clear all clc I=0; %赋予I20的初始值 for n=0:19; I=(-1/5)*I+1/(5*(20-n)); %有递推公式得 end fprintf('I0=%f\n',I)M=10000; for i=0:19; M=(-1/5)*M+1/(5*(20-i));%有递推公式得 end fprintf('M0=%f\n',M) 输出结果为： I0=0.182322 M0=0.182322（3）由输出结果可看出第一种算法为不稳定算法，第二中算法为稳定算法。
由于误差*000***21111120115(5)5()555nn n n n n n n n n e I I e I I I I I I e e e n n------=-=-=-+--+=-===第一种算法为正向迭代算法，每计算一步误差增长5倍，虽然所给的初始值很接近，随着n 的增大，误差也越来越大。
数值分析计算实习第一题
直接用定义: ������������(������������)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
数值计算方法第一次上机实习报告
数值计算方法第一次实习报告一、 实习题目。
1. 分别用高斯-赛德尔迭代法、雅克比迭代法、列主元消去法解下列方程组。
⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++ 1.84711.2671x 0.2568x 0.2471x 0.2368x1.74710.2271x 1.2168x 0.2071x 0.1968x 1.64710.1871x 0.1768x 1.1675x 0.1582x 1.54710.1490x 0.1397x 0.1254x 1.1161x4321432143214321 2. 用迭代法求x 5-x-0.2=0的正根，要求精确到小数点后第五位。
二、算法原理。
1、雅可比迭代法基本原理 将矩阵分解为，其中则式可记为，变形可得，可逆时，有于是得到迭代的过程为式中，，即2、高斯-赛德尔迭代法基本原理赛德尔迭代法是对雅可比迭代法的一种改进，雅可比迭代法是在每一步计算的各个分量时均只用到中的分量。
实际上，在计算时，分量都已经计算出来而没有被直接利用，因此可以考虑以来代替计算。
即121311121232222121-1,212121000, ,0000n n n n nn a a a a a a a a D L U aa a a a a a ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==-=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦=Ax b ()=D L U x b--()=+Dx L U x b+D ()11220nn a a a ≠ ()11=D +x L U x D b --+=+x Bx f ()11D,=D B L U f b --=+()=11=+=1,2,,n i i ij j j ii j i x b a x i n a ≠⎛⎫ ⎪ ⎪ ⎪⎝⎭∑ ()+1k x ()k x ()+1k i x ()()+1+1-1,,k k i i x x ()()+1+1-1,,k k i i x x ()()1-1,,kki x x矩阵形式为，可得，于是赛德尔迭代法的矩阵形式为式中，。
[整理]-8-《数值计算方法》实习作业(模板-小)
![[整理]-8-《数值计算方法》实习作业(模板-小)](https://img.taocdn.com/s3/m/dfc299b7d4d8d15abe234ef2.png)
2.1函数图形与极限2.1.1 实验目的1.熟悉Mathematica 基本绘图语句。
2.掌握函数极限的有关操作命令。
3.学会利用Mathematica 软件对函数进行分析研究。
4.熟悉Mathematica 二元函数绘图语句。
2.1.2 实验内容【基本语句】1.Plot[f[x],{x,xmin,xmax},选项]; 功能: 画出函数f[x] 从min 到max 间的图形;2.Plot[{f1[x],f2[x],...},{x,xmin,xmax},选项]; 功能: 在同一坐标系下画出函数f1,f2,...的图形。
3. ParametricPlot[{fx,fy},{t,tmin,tmax}]; 功能: 画出参数方程fx=x(t),fy=y(t)的图形;ParametricPlot[{{f1x,f1y},{f2x,f2y}},{t,tmin,tmax}]; 功能：在同一坐标系下画出用参数方程表示的两幅函数图形。
【备注】fx,fy 的给出方式:⑴fx=x(t) , fy=y(t)⑵fx=x ,fy=f(x)与fx=f(x) ,fy=x 构成反函数的图形关系⑶r=r(t) , fx=r(t)Cos(t) , fy=r(t)Sin(t)4. Show[tu1,tu2] 功能：将tu1及tu2两幅函数图形重叠在一起，将两个函数图形一起显示。
5. Plot3D[f[x,y],{x,x0,x1},{y,y0,y1}] 功能：作出函数f[x,y]在区域[x0,x1]×[y0,y1]上的图形； ParametricPlot3D[{x[u,v],y[u,v],z[u,v]},{u,u0,u1},{v,v0,v1}] 功能：作出参数方程表示的曲面。
6. Limit[f[x],x->x0] 功能：求函数f[x]在x0处的极限。
7. Limit[f[x],x->x0,Direction->+1] 功能：求函数f[x]在x0处的左极限。
数值计算方法上机实验题
数值计算方法上机实验实验内容：1.要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.2．给定积分dx e x⎰31和dx x ⎰311 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1）变步长梯形法; 2）变步长 Simpson 法; 3) Romberg 方法.算法描述：1、复合梯形法：⎰=tdt t a t V 0)()( ))()(2)((211∑-=++=n k k n b f x f a f hT输入 被积函数数据点t,a. 输出 积分值.n T复合Simpson 法：⎰=tdt t a t V 0)()( ))()(2)(4)((6101121∑∑---=++++=n k n k k k n b f x f x f a f hS输入 被积函数f(x),积分区间[a,b]和n 输出 复合Simpson 积分值n S步1 .);()(;a x b f a f S nab h n ⇐-⇐-⇐ 步2 对n k ,,2,1 =执行).(2;2);(4;2x f S S hx x x f S S h x x n n n n +⇐+⇐+⇐+⇐步3 n n S hS ⨯⇐6步4 输出n SRomberg 积分法：根据已知数据对其进行多项式拟合得出p(x);f(x)⇐p(x); 输入 被积函数f(x)，积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111⇐===+⇐-⇐k R C S b f a f hT a b h 步2 反复执行步3→步9. 步3 .2;0h a x S +⇐⇐ 步4 反复执行步5→步6. 步5 ;);(h x x x f S S +⇐+⇐步6 若x ≥b,则退出本层循环. 步7 执行.6316364;1511516;3134;2212212212212C C R S S C T T S S h T T -⇐-⇐-⇐+⇐步8 执行.1;;;;;2;2121212112+⇐⇐⇐⇐⇐⇐-⇐k k R R C C S S T T hh R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ⇐ 步11 输出.2n R2、变步长梯形算法：功能 求积分⎰ba)(dx x f ，允许误差为ε。
数值计算方法上机题
习题二问题：1．编制通用子程序对n+1个节点x i及y i=f(x i) (i=1,…n)（1）n次拉格朗日插值计算公式；（2）n次牛顿向前插值计算公式；（3）n次牛顿向后插值计算公式；（一）程序流程图（1）拉格朗日插值程序流程图（2）牛顿向前插值程序流程图（3）牛顿向后插值程序流程图。
（二）源程序见主程序清单问题：2．计算（1）已知f(x)=lnx,，[a，b]=[1，2]，取h=，x i=1+ih，i=0，1， (10)用通用程序（1），（3）计算及的近似值；（一）程序清单/* program of question , page 61 */#include ""#include ""main(){ int i,flag=0;double z1,z2,x[11],y[11],t,s1,s2,z[10],c[11][11],log(double),ntb(),L(); for(i=0;i<=10;i++){x[i]=1+*i;y[i]=log(x[i]);}printf("data x:\n");for(i=0;i<=10;i++){flag++;printf("%11.6f",x[i]);if(flag%4==0)printf("\n");}printf("\ndata y:\n");flag=0;for(i=0;i<=10;i++){flag++;printf("%11.6f",y[i]);if(flag%4==0)printf("\n");}printf("\nThe true value:\n");printf(" =%f =%f\n",log,log);z1=L(x,y,10,;z2=L(x,y,10,;t=[10])/;s1=ntb(y,10,t,z,c);s2=ntb(y,10,t,z,c);t=[10])/;printf("The approximate value:\n");printf(" L=%f L=%f\n",z1,z2);printf(" NTB=%f NTB=%f\n",s1,s2);}double L(double x[],double y[],int n,double t){ int i,k; double z=,s;if(n==1)z=y[0];for(k=0;k<=n;k++){ s=;for(i=0;i<=n;i++) if(i!=k)s=s*(t-x[i])/(x[k]-x[i]);z=z+s*y[k]; }return z;}double ntb(double y[],int n,double t,double z[],double c[][11]) { int i,j,sn=n;double s;z[0]=t;for(i=1;i<=n-1;i++) z[i]=z[i-1]*(t+i)/(i+1);for(i=0;i<=n;i++) c[i][0]=y[sn--];for(j=1;j<=n;j++)for(i=0;i<=n-j;i++) c[i][j]=c[i][j-1]-c[i+1][j-1];s=y[n];for(i=0;i<=n-1;i++) s=s+z[i]*c[0][i+1];return s;}(二)运行结果data x:data y:The true value:= =The approximate value:L=0.431782 L=NTB= NTB=问题：（2）f(x)=1/(1+25x2), |x|≤1取等距节点n=5和n=10，用通用程序（1），（2）依次计算x=+ih（i=0，1，…，19，h=）处f(x)的近似值，并将其结果与其真实值相比较。
数值计算方法上机实习题答案.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） 分析绝对误差。
数值计算方法上机实习题NEW
数值计算方法上机实习题1． 设，（1） 由递推公式，从 , 出发，计算 ； 程序如下：function I=myhs(I0,n) if n>=1I=myhs(I0,n-1)*(-5)+1/n; elseif n==0 I=I0; end命令行窗口输入： I0=0.1822;I1=myhs(I0,20); I0=0.1823;I2=myhs(I0,20);运行结果：当I0=0.1822时，I20=-1.1593e+10。
当I0=0.1823时，I20= -2.0558e+9。
（2） ，20=10000I , 用，计算0I ； 程序如下：function I=myhs2(I20,n) if (n<20)I=myhs2(I20,n+1)*(-1)/5+1/(5*(n+1)); elseif n==20 I=I20; end命令行窗口输入：I20=0;I1=myhs2(I20,20); I20=10000;I2=myhs2(I20,20);运行结果：当I 20=0时，I 0=0.182321556793955。
当I 20=10000时，I 0= 0.182321556898812。
（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
根据上式，假设I n 的真值为I*，误差为e n ，即e=I*-I n 。
综合递推式，有e n =-5*e n-1，这意味着哪怕开始只有一点点误差，只要n 足够大，按照这种计算一步误差增长5倍的方式，所得的结果总是不可信的，因此整个算法是不稳定的。
根据上式，假设I n 的真值为I*，误差为e n ，即e=I*-I n 。
综合递推式，有e n-1. =(-1/5)*e n ，按照这种计算误差会以每步缩小到1/5的方式进行，根据（2）得到的结果而言，该算法是相对稳定的。
2． 求方程0210=-+x e x 的近似根，要求41105-+⨯<-k k x x ，并比较计算量。
数值方法计算实习题1
信计091 龚立丽200900901004数值方法计算实习题要求：1、用Matlab语言或你熟悉的其他算法语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。
一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和20次Lagrange插值多项式对数据进行插值。
用图示出给定的数据，以及()s x和20()L x。
x0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0y 1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25x7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3y 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25解：>> x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 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）三次样条插值法在MATLAB中编写m文件function[f,f0]=scyt(x,y,y2_1,y2_N,x0)%y2_1和y2_N分别为自然边界条件%插值点x的坐标：x0syms t;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等');return;endfor i=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);for i=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i) );A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=zgf(A,c);h=x(index+1)-x(index);f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1 )*(x(index+1)-t)*(t-x(index))^2/h/h;f0=subs(f,'t',x0);其中的zgf函数为追赶法，其程序为function x=zgf(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0！');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1); end在MATLAB 中输入指令 >> [f,f0]=scyt(x,y,0,0) f =1000/729*(27/5*t-1377/100)*(t-39/10)^2+1000/729*(522/25-24/5*t)*(t-3)^2+100/81*(-6396162352027119/288230376151711744*t+19188487056081357/288230376151711744)*(39/10-t)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*t)*(t-3)^2 f0 =2.5851 得三次样条插值函数 S(x)=1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(x-3)^2+100/81*(-6396162352027119/288230376151711744*x+19188487056081357/288230376151711744)*(39/10-x)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*x)*(x-3)^2>> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> title('试验一--三次样条插值图示')024********0.511.522.53试验一--三次样条插值图示（2）用拉格朗日法插值 %定义Lagrange 程序function f=Language(x,y,x0) syms t ;if (length(x)==length(y)) n=length(x);elsedisp('x 和y 的维数不相等'); return ; end f=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 ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦，则方程组A x b =有准确解[]1111Tx =。
数值计算方法上机实验题
数值计算方法上机实验实验内容：1.要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.2．给定积分dx e x⎰31和dx x ⎰311 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1）变步长梯形法; 2）变步长 Simpson 法; 3) Romberg 方法.算法描述：1、复合梯形法：⎰=tdt t a t V 0)()( ))()(2)((211∑-=++=n k k n b f x f a f hT输入 被积函数数据点t,a. 输出 积分值.n T复合Simpson 法：⎰=tdt t a t V 0)()( ))()(2)(4)((6101121∑∑---=++++=n k n k k k n b f x f x f a f hS输入 被积函数f(x),积分区间[a,b]和n 输出 复合Simpson 积分值n S步1 .);()(;a x b f a f S nab h n ⇐-⇐-⇐ 步2 对n k ,,2,1 =执行).(2;2);(4;2x f S S hx x x f S S h x x n n n n +⇐+⇐+⇐+⇐步3 n n S hS ⨯⇐6步4 输出n SRomberg 积分法：根据已知数据对其进行多项式拟合得出p(x);f(x)⇐p(x); 输入 被积函数f(x)，积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111⇐===+⇐-⇐k R C S b f a f hT a b h 步2 反复执行步3→步9. 步3 .2;0h a x S +⇐⇐ 步4 反复执行步5→步6. 步5 ;);(h x x x f S S +⇐+⇐步6 若x ≥b,则退出本层循环. 步7 执行.6316364;1511516;3134;2212212212212C C R S S C T T S S h T T -⇐-⇐-⇐+⇐步8 执行.1;;;;;2;2121212112+⇐⇐⇐⇐⇐⇐-⇐k k R R C C S S T T hh R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ⇐ 步11 输出.2n R2、变步长梯形算法：功能 求积分⎰ba)(dx x f ，允许误差为ε。
数值计算方法上机实习题考证
---------------------------------------------------此文档包含我们计算方法的经典算法包含（数值计算方法上机实习题）1． 设⎰+=105dx x x I nn ， （1） 由递推公式nI I n n 151+-=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用n I I n n 51511+-=-，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
(1) 解答：n=0，0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=⎰⎰⎰x d x dx x dx x x I nn 这里可以用for 循环，while 循环，根据个人喜好与习惯：for 循环程序： While 循环程序：I=0.1823; I=0.1823;for n=1:20 i=1;I=(-5)*I+1/n; while i<21End I=(-5)*I+1/i;I i=i+1;fprintf('I20=%f',I) endI = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009(2) 粗略估计I 20： Mathcad 计算结果： for 循环程序： While 循环程序： >> I=0.007998; I=0.007998;>> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n);>> I n=n+1;I =0.0083 end>> II =0.0083(3) 算法误差分析：计算在递推过程中传递截断误差和舍入误差第一种算法：（从1——>20）01x x 205x +⎛⎜⎜⎜⎠d 7.998103-⨯=*000e I I =-***21111120115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===误差放大了5n 倍，算法稳定性很不好；第二种算法：（从20——>1）*n n ne I I =- ***111111111()()555555n n n n n n n n e I I I I I I e n n ---=-=-+--+=-=0111...()55n n e e e === 误差在逐步缩小，算法趋近稳定，收敛。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算方法上机实习题
1． 设⎰+=1
05dx x
x I n
n ， （1） 由递推公式n
I I n n 1
51+
-=-，从I 0=0.1824, 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n
I I n n 51
5111+-
=--，计算0I ；
（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
答：第一个算法可得出
e 0=|I 0−I 0
∗| e n =|I n −I n ∗|=5n |e 0|
易知第一个算法每一步计算都把误差放大了5倍，n 次计算后更是放大了5n 倍，可靠性低。
第二个算法可得出
e n =|I n −I n ∗| e 0=(15
)n
|e n |
可以看出第二个算法每一步计算就把误差缩小5倍，n 次后缩小了5n 倍，可靠性高。
2． 求方程0210=-+x e x
的近似根，要求41105-+⨯<-k k x x ，并比较计算量。
（1） 在[0，1]上用二分法；
计算根与步数程序：
fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x;
f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1);
fprintf('root=%6.8f ,n=%d \n',root,n);
计算结果显示：
root=0.09057617 ,n=11
（2） 取初值00=x ，并用迭代10
21
x k e x -=+；
（3） 加速迭代的结果;
（4） 取初值00 x ，并用牛顿迭代法；
（5）分析绝对误差。
答：可以看到，在同一精度下，二分法运算了11次，题设迭代算式下运算了4次，加速迭代下运算了2次，牛顿迭代下运算了2次。
因不动点迭代法和二分法都是线性收敛的，但二分法压缩系数比题设迭代方法大，收敛速度较慢。
加速迭代速度是超线性收敛，牛顿法是二阶，收敛速度快。
3．钢水包使用次数多以后，钢包的容积增大，数据如下：
试从中找出使用次数和容积之间的关系，计算均方差。
（用
ax b
y
c x
+
=
+
拟合）
构造函数子程序：
function delta=delta1(f,x,y)
a=f(1);
b=f(2);
c=f(3);
delta=0;
for k=1:15
delta=delta+((x(k)+c)*y(k)-(a*x(k)+b))^2;
end
计算结果显示：
拟合出的方程为:(x+-0.7110)y=11.2657x+-15.5024
均方差为:0.33165089
总结：
指标选择，因题设方程为非线性的，要转化为线性方程故需提指标为：
其驻点方程为:计算结果显示：
4．设⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫
⎝⎛----------------=410100141010014101101410010141001014A ，⎪⎪⎪⎪⎪⎪
⎪⎪⎭
⎫
⎝⎛--=625250b ，b x =A
分析下列迭代法的收敛性，并求42
110-+≤-k k x x 的近似解及相应的迭代次数。
（2） GAUSS-SEIDEL 迭代；
5．用逆幂迭代法求⎪⎪⎪⎭
⎫ ⎝⎛=111123136A 最接近于11的特征值和特征向量，准确到3
10-。
6．用经典R-K 方法求解初值问题
（1）⎩⎨⎧-+-='++-='x x y y y x y y y sin 2cos 22sin 22212211
，]10,0[∈x ，
⎩⎨
⎧==3
)0(2
)0(21y y ；
（2）⎩⎨⎧-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212
211
，]10,0[∈x ，
⎩⎨
⎧==3)0(2
)0(2
1y y 。
和精确解⎩⎨⎧+=+=--x
e x y x
e x y x
x cos 2)(sin 2)(21比较，进行误差分析得到结论，图形显示精确解和数值解。
计算结果显示：
(1)Ode 23的题1
(2)Ode 3的题2
(3)Ode 45的题1
(4)Ode 45的题2
7．用有限差分法求解边值问题（h=0.1），并图形显示。
⎩⎨
⎧==-=+-''1
)1()1(0
)1(2y y y x y .。