计算数值方法上机作业

合集下载

数值计算方法上机实习题考证

数值计算⽅法上机实习题考证--------------------------------------------------- 此⽂档包含我们计算⽅法的经典算法包含（数值计算⽅法上机实习题）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) ................................................................................................ 错误!未定义书签。

数值计算方法上机题目资料

，在程序设计时要充分考虑哪些变量应该可变的。
2. 程序输入、输出用文件形式。 3. 编程语言要求用C, 编程环境TC或VC 4. 程序要求调试通过。 5. 每个方法要求给出一个具体的算例（可选
作业题）来验证。
五、上机报告要求
1.报告内容包括：
每种方法的算法原理及程序框图。程序使用说明具体算例及结果
上机调试体会及收获。
2.报告要手写。
六、上机报告及源程序提交时间
1.上机报告在考试当天提交。 2.源程序在考试前提交。
提交格式：文件夹（班级+姓名）
输入文件程序文件夹输出文件
不要拷贝其它文件！！！源程序
六、上机报告及源程序提交时间
源程序提交：把以上文件压缩后，发送到以下邮箱：
haoyq@
七、考核方式
1.算法手算笔试（80%）+上机内容笔试（10%）+上机报告（10%）
2.上机内容笔试可能形式：
编一段算法程序给出一段算法程序，说明算法的名称。程序填空程序改错（包括算法和语法的错误）
数值计算方法上机练习
一、上机练习目的
复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
二、上机练习任务
• 利用计算机基本C语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
• 掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
• 写出上机练习报告。
三、数值计算方法上机题目

数值代数与计算方法上机作业

数值代数与计算方法上机作业作业一：Matlab的基本操作P311.根据习题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的解法P401. 使用程序2.1求解下面每个函数的不动点（尽可能多）近似值，答案精确到小数点后12为。

同时，构造每个函数和直线y=x来显示所有不动点。

（a）（b）（c）（d）g(x)?x5?3x3?2x2?2g(x)?cos(sin(x)) g(x)?x2?in(x?0.15) g(x)?xx?cos(x)P493. 修改程序2.2和程序2.3，使得输出分别类似于表2.1和表2.2的矩阵（即矩阵的第一行应当为[0a0 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.作业四：插值与多项式逼近P154Matlab的矩阵特性使其能够快速计算一个函数在其多个点处的值。

例如，如果X=[-1 0 1],则sin(X)将得到[sin(-1),sin(0),sin(1)]。

数值计算方法上机作业

《数值计算方法》上机作业学院：机械学院专业：机械设计及理论姓名：陈国保学号：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）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；（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.二分法2.不动点迭代法三、结构程序设计代码1.二分法1).定义所求解函数function [ y ] = f( x )y = x^3 + 4*x^2 - 10;end2).执行算法%初始化，设置区间端点a、b，误差限tola = 1;b = 2; tol = 0.5*10^(-6);k = 0; fa = f(a);%设置最大二分次数为30for k = 1:50p = (a + b)/2; fp = f(p);if(fp == 0 || (b - a)/2 < tol)breakendif(fa * fp < 0)b = p;elsea = p;enddisp('近似解p = ');disp(vpa(p,10)); disp('迭代次数k = ');disp(k);end2.不动点迭代法1).定义不动点方程g(x)function [ y ] = g( x )y = x^3 + 4*x^2 + x - 10;end2).执行算法%初始化，设置误差限，设置初值p0tol = 0.5*10^(-6);k = 1; p0 = 1.5;%迭代次数为10次while k <= 10p = g(p0);if abs(p - p0) < tolbreakenddisp('近似解p = ');disp(vpa(p,10)); disp('迭代次数k = ');disp(k);k = k + 1; p0 = p;end四、结果及其讨论1.二分法结果由于结果较长，只取了一部分，从图中可以看出，迭代20次可得到误差限范围内的近似解p=1.36522。

数值计算方法上机1

插值1. 用系数矩阵求方程组的方法求得拉氏插值多项式的方法在MATLAB 中，输入命令format rat %使数据格式采用有理数形式x=[-1,1,2,5]’,y=[-7,7,-4,35]’ %输入列向量A=[x.^0,x.^1,x.^2,x.^3] %构造出线性方程组(2.2)的系数矩阵A\yans=[10,5,-10,2]’(运行文件名统一取为ff1.m),其结果表明 23()105102p x x x x =+-+2．拉格朗日插值函数的计算机实现编制计算拉格朗日插值的M 文件：以下是拉格朗日插值的名为y_lagrl 的M 文件：function y=y_lagr1(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end运用此文件到具体插值问题中的例子：选择函数y=exp(-x 2) (-2≤x ≤2),在n 个节点上（n 不要太大，如5~11）用拉格朗日插值方法，计算m 个插值点的函数值（m 要适中，如50~100），并画出插值图形。

n=7;m=61;x=-2:4/(m-1):2;y=exp(-x.^2);z=0*x;x0=-2:4/(n-1):2;y0=exp(-x0.^2);y1=y_lagr1(x0,y0,x);plot(x,y,’r ’,x,y1,'b:')；（文件名取为ff2.m ）运行后结果：3．龙格现象在MA TLAB 中，通过编程可以观察到这种“龙格现象”。

%步1 定义名为f.m 的函数文件function y=f(x)y=1./(1+x.^2);%步2 编写计算拉格朗日插值多项式值并作图程序%作被插函数图象x=-5:1:5;y=f(x);plot(x,y);%固化图形屏幕Hold on%计算下列节点处插值多项式的值z=-5:0.005:5;m=max(size(z));n=10;for k=1:mL(k)=0;for i=1:n+1t(i)=1;for j=1:n+1if j~=it(i)=t(i)*(z(k)-x(j))/(x(i)-x(j));endendL(k)=L(k)+t(i)*y(i);endendplot(z,L,'b')hold off在matlab 中实现，（文件名为ff3.m ）结果为：4．分段线性插值函数逼近被插函数的计算机摸拟作函数211xy +=在节点),,1,0(105n i i n x i =⋅+-=的分段线性插值图象，并将它们进行比较，观察可发现，随着n 的增大，两者吻合得越来越好，龙格现象并未发生。

数值计算方法上机答案

2
0.23

22.322

b

54 240.236

29.304

117.818
2.991 1 13
3.907 0.017
3
1.007 2.1 0.5 61.705 0.91 4.5
1.006 6.3 1
12.17 4.918
[ 1.0, 2.0, 3.0, 4.5, 5.0, 21.8]
ans =
1.0000 -2.0000
2.9999 -4.0001
5.0000 -6.0008
17 1 4 3 1 2 3 7

2
10 1
7
2 1
1
4
1 1 8 2 5 2 1 1
3、
A

max1=0; for i=k:n
if abs(A(i,k))>max1 max1=abs(A(i,k));r=i;
end end if max1<1e-6
flag='failure';return; end if r>k
for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
[ -1.0, 34.211, -1.0, -2.1, 6.3, -1.7]
[
0, 0.5, 13.0, -0.5, 1.0, -1.5]
[ 4.501, 3.11, -3.907, -61.705, 12.17, 8.999]
[ 0.101, -0.812, -0.017, -0.91, 4.918, 0.1]
k x(k) 0 1.500000 1 0.784472 2 0.739519

数值计算方法上机实验题

数值计算方法上机实验实验内容：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)的近似值，并将其结果与其真实值相比较。

计算方法上机题目

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

数值计算方法上机实习题答案.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-1解析

数值计算方法上机题目11、实验1. 病态问题实验目的：算法有“优”与“劣”之分，问题也有“好”和“坏”之别。

所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。

希望读者通过本实验对此有一个初步的体会。

数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。

病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。

问题提出：考虑一个高次的代数多项式∏=-=---=201)()20)...(2)(1()(k k x x x x x p （E1-1）显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。

现考虑该多项式方程的一个扰动0)(19=+xx p ε (E1-2)其中ε是一个非常小的数。

这相当于是对（E1-1）中19x 的系数作一个小的扰动。

我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。

实验内容：为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数u ＝roots （a ）其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。

设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程0...1121=++++-n n n n a x a x a x a的全部根，而函数b=poly(v)的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。

可见“roots ”和“Poly ”是两个互逆的运算函数.ve=zeros(1,21); ve(2)=ess;roots(poly(1:20))+ve)上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。

实验要求：（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。

数值分析上机作业(2)

一、数值求解如下正方形域上的Poisson 方程边值问题 2222(,)1,0,1(0,)(1,)(1),01(,0)(,1)0,01u u f x y x y x y u y u y y y y u x u x x ⎧⎛⎫∂∂-+==<<⎪ ⎪∂∂⎪⎝⎭⎨==-≤≤⎪⎪==≤≤⎩二、用椭圆型第一边值问题的五点差分格式得到线性方程组为2,1,1,,1,10,1,,0,141,?,?,?,?0,1i j i j i j i j i j ijj N j i i N u u u u u h f i j N u u u u i j N -+-+++----=≤≤====≤≤+，写成矩阵形式Au=f 。

其中1．三、编写求解线性方程组Au=f 的算法程序, 用下列方法编程计算, 并比较计算速度。

2．用Jacobi 迭代法求解线性方程组Au=f 。

3．用块Jacobi 迭代法求解线性方程组Au=f 。

4．用SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

1122N N v b v b u f v b ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭4114114ii A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭11,12,1,121,22,2,21,2,,2211,12,1,121,22,2,221,2,,(,,...,),(,,...,),......,(,,...,)(,,...,)?,(,,...,)?,......,(,,...,)?1,999,0.10.011T T N N TN N N N N T T N N T N N N N N v u u u v u u u v u u u b h f f f b h f f f b h f f f h N h N ====+=+=+===+取或则或,1,,1,2,...,i j f i j N== 1122NN A I I A A I I A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭5．用块SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

数值计算方法上机实验题

数值方法上机作业

第二章插值法1.通用程序编制（1）n次拉格朗日插值计算公式#include<stdio.h>void main(){int n,m,i,j,k;float a[100],b[100],c,y[100],x[100]; /*输入插值多项式次数*/scanf("%d",&n);/*输入所求未知量个数*/scanf("%d",&m);/*输入插值节点坐标*/for(i=0;i<=n;i++){scanf("%f",&a[i]);scanf("%f",&b[i]);}/*输入所求点x坐标*/for(i=0;i<m;i++){scanf("%f",&x[i]);}/*计算*/for(j=0;j<m;j++){y[j]=0;for(k=0;k<=n;k++){c=1;for(i=0;i<=n;i++){if(k!=i)c=c*(x[j]-a[i])/(a[k]-a[i]);}c=c*b[k];y[j]=y[j]+c;}}/*输出计算结果*/for(i=0;i<m;i++){printf("%f\n",y[i]);}getchar();}（2）n次牛顿向前插值计算公式#include<stdio.h>void main(){int n,m,i,j;float x,h,a,b[100],c[100],y[100],t[100]; /*输入插值多项式次数*/scanf("%d",&n);/*输入所求未知量个数*/scanf("%d",&m);/*输入步长*/scanf("%f",&h);/*输入xo点坐标*/scanf("%f",&x);/*输入插值节点坐标*/for(i=0;i<=n;i++){scanf("%f",&b[i]);}/*输入所求点t值*/for(i=0;i<m;i++){scanf("%f",&t[i]);}/*造向前差分表*/for(i=0;i<n;i++){for(j=i;j<n;j++){c[j+1]=b[j+1]-b[j];}for(j=i+1;j<=n;j++){b[j]=c[j];}}/*计算*/for(i=0;i<m;i++){a=1;y[i]=0;for(j=0;j<=n;j++){y[i]=y[i]+a*b[j];a=a*(t[i]-j)/(j+1);}}/*输出计算结果*/for(i=0;i<m;i++){printf("x=%f y=%f\n",x+t[i]*h,y[i]); }getchar();}（3）n次牛顿向后插值公式#include<stdio.h>void main(){int n,m,i,j;float x,h,a,b[100],c[100],y[100],t[100];/*输入插值多项式次数*/scanf("%d",&n);/*输入所求未知量个数*/scanf("%d",&m);/*输入步长*/scanf("%f",&h);/*输入Xn点坐标*/scanf("%f",&x);/*倒着输入插值节点坐标*/for(i=0;i<=n;i++){scanf("%f",&b[i]);}/*输入所求点t值*/for(i=0;i<m;i++){scanf("%f",&t[i]);}/*造向后差分表*/for(i=0;i<n;i++){for(j=i;j<n;j++){c[j+1]=b[j]-b[j+1];}for(j=i+1;j<=n;j++){b[j]=c[j];}}/*计算*/for(i=0;i<m;i++){a=1;y[i]=0;for(j=0;j<=n;j++){y[i]=y[i]+a*b[j];a=a*(t[i]+j)/(j+1);}}/*输出计算结果*/for(i=0;i<m;i++){printf("x=%f y=%f\n",x+t[i]*h,y[i]); }getchar();}2.计算（1）运用通用程序1计算结果如下x=1.54时, y=0.431808x=1.98时, y=0.683157运用通用程序3计算结果如下x=1.54时, y=0.431808x=1.98时, y=0.683158 （2）n=5时，运用通用程序1计算结果如下运用通用程序2计算结果如下n=10时，运用通用程序1计算结果如下运用通用程序2计算结果如下计算结果与真实值的比较与真实值相比，在[-0.15 0.15]的范围内，n=10比n=5能更好的接近真实值。

数值计算方法上机实验报告1

数值计算方法上机实验报告1华北电力大学上机实验报告课程名称：数值计算方法专业班级：学生姓名：学号：指导教师：张建成实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。

上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。

掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。

一、列主元素消去法求解线性方程组 1、算法原理为避免绝对值很小的元素作为主元，在每次消元之前增加一个选主元的过程，将绝对值大的元素交换到主对角线的位置。

列主元素消元法是当变换到第k 步时，从k 列的kk a 及以下的各元素中选取绝对值最大的元素，然后通过二交换将其交换到kk a 的位置上。

2、输入输出变量ija ：为系数矩阵的各个系数K ：表示到第k 步消元 3、具体算例输入增广矩阵为： 3 1 2 -3 8 2 1 3 22 3 2 1 28解得：1x =6，2x =4，3x =2；二、LU 分解法求解线性方程组1、算法原理应用高斯消去法解n 阶线性方程Ax b =经过1n -步消去后得出一个等价的上三角形方程组()()n n A x b =，对上三角形方程组用逐步回代就可以求出解来。

这个过程也可通过矩阵分解来实现。

将非奇异阵分解成一个下三角阵L 和上三角阵U 的乘积A LU =称为对矩阵A 的三角分解，又称LU 分解。

根据LU 分解,将Ax b =分解为Ly bUx y =??=?形式,简化了求解问题。

2、输入输出变量ij a 为系数矩阵元素i b 为常数矩阵系数,i j i jl u 分别为下、上三角矩阵元素 k 表示第k 步消元 3、具体算例输入增广矩阵 3 2 3 4 39 3 -2 2 14 4 2 3 43 解得： 6 5 3三、拉格朗日插值1、算法原理设函数()y f x =在区间[a,b]上有节点01,,,,n x x x 上的函数值，构造一个次数不超过n次的代数多项式1110()n n n n p x a x a x a x a --=++++ ,使 (),0,1,,i i P x y i n == 。

数值计算方法上机实习题

数值计算方法上机实习题(总13页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--数值计算方法上机实习题1．设⎰+=105dx xx I nn ，（1）由递推公式nI I n n 151+-=-，从I 0=0.1824, 0=0.1823I 出发，计算20I ；（2） 20=0I ，20=10000I , 用nI I n n 515111+-=--，计算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=(1)n|e n |可以看出第二个算法每一步计算就把误差缩小5倍，n 次后缩小了5n 倍，可靠性高。

2．求方程0210=-+x e x 的近似根，要求41105-+⨯<-k k x x ，并比较计算量。

（1）在[0，1]上用二分法； (1) [0，1]上的二分法二分法子程序：function [root,n]=EFF3(f,x1,x2) %第二题(1)二分法f1=subs(f,symvar(f),x1);%函数在x=x1的值 f2=subs(f,symvar(f),x2);%x=x2 n=0;%步数er=5*10^-4;%误差 if(f1==0) root=x1; return; elseif(f2==0) root=x2; return;elseif(f1*f2>0)disp('两端点函数值乘积大于0！'); return; elsewhile(abs(x1-x2)>er)%循环 x3=(x1+x2)/2;f3=subs(f,symvar(f),x3); n=n+1; if(f3==0) root=x3; break;elseif(f1*f3>0) x1=x3; elsex2=x3; end endroot=(x1+x2)/2;%while 循环少一步需加上 end计算根与步数程序：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=% ,n=%d \n',root,n);计算结果显示：root= ,n=11（2）取初值00=x ，并用迭代1021x k e x -=+；(2) 初值x 0=0迭代迭代法子程序：function [root,n]=DDF(g,x0,err,max) （接下页） %root 根，n+1步数，g 函数，x0初值，err 误差，max 最大迭代次数 X(1)=x0; for n=2:maxX(n)=subs(g,symvar(g),X(n-1)); c=abs(X(n)-X(n-1)); root=X(n); if(c<err)计算根与步数程序：syms x;f=(2-exp(x))/10; （接下页） x0=0;err=5*10^(-4); max=100;[root,n]=DDF(f,x0,err,max); fprintf('root=% ,n=%d \n',root,n);计算结果显示： root= ,n=4（3）加速迭代的结果;（4）取初值00 x ，并用牛顿迭代法；（5）分析绝对误差。