计算方法上机题目
目录
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 u u
=-,使得H x
v
σ=，其中σ为常数，当取单位向量2
x v u
x v
σσ-=
-时，由u 确定的矩阵H 必定满足H x
v
σ=，所以在计算过程中取u 的值为上述值。
设A 是一个()
m n m
n ?≥阶矩阵且它的列向量线性无关，则利用豪斯霍尔德
变换可以把A 逐步化为上梯形矩阵，设
()11
12121
22
212
1
2
,,
,n n
n
m m m n a a a a a a A a a a a a a ??
?
?== ? ???
具体变换过程如下：
设()1,2,
,i e i
n =是
m 维单位坐标向量。
第1步 为把矩阵A 的第一列
()
01
a 化为()1,
0,
,0T
σ，
取
()
()01
1,1,0,,0T
m
x a v e R
===∈ （或取1v
e =-），根据上式可得，取
()
()
01111101
2
1
11
2
a e u a e σωωσ-=
=
-
其中
()
()
001
21
11
21111, a e a e ωσωσ=-=-
()
(
)
()
(
)011
111111111
2
01
2
1
1111
22, T T
T
T
u u
a a ωωωωωωασσωασσ=
=
=
=--
令
1
111
111
22T
T
m m H I u u I αωω-=-=-，用1H 左乘()0A 得，
()
()
()
()
()
()
()
()
()
(
)
1000011
1
12
11111123
,,
, =,,,
,n n
A
H A
H a H a H a e a a a σ==
程序框图
2
2i i i i i i
T
a e I u u
σ--
计算实习
用豪斯霍尔德变换求方程组x
A b = ，其中
54756753941287886537810987756, 5
79119755368891089587877810105756
7
5
9
10
1052A b ????????????
??
??????==????
??
??????????
?????
???
Matlab 代码 %使用说明：
%需自己输入矩阵A 及b 的值
%变量Q,R 分别为进行QR 分解后的结果 clear clc format long load('A 矩阵.mat') load('b 矩阵.mat')
%调用函数qrhs 进行QR 分解 [Q,R]=qrhs(A); [~,n]=size(R);
fprintf('您输入的矩阵阶数'); n y=Q'*b;
%回代过程
x(n,1)=y(n,1)/R(n,n); for i=n-1:-1:1 s=y(i,1);
for j=i+1:n
s=s-R(i,j)*x(j,1);
end
x(i,1)=s/R(i,i);
end
x
被调用函数qrhs
function[Q,R]=qrhs(A)
format long
[~,n]=size(A);
Q=eye(n);
for j=1:n-1
B=norm(A(j:n,j));
Y=zeros(n-j+1,1);
Y(1,1)=-sign(A(1,j))*B;
X=A(j:n,j);
I=eye(n-j+1);
N=I-(2/(norm(X-Y))^2)*(X-Y)*(X-Y)';
H=[eye(j-1) zeros(j-1,n-j+1);zeros(n-j+1,j-1) N];
A=H*A;
Q=Q*H;
end
R=A;
Q;
R;
end
3. 共轭梯度法求解线性方程组
计算原理
当A 是n 阶对称正定矩阵，若*x 是二次函数()12
T T
f x x A x b x
=-的极小值点，
则*x 是方程组A x
b =的解，即
()()n
*
*
x R
m in A x b f
x f
x ∈=?=
共轭梯度法在形式上具有迭代法的特征，即给定初始向量(0)x ，由迭代格式
()
()
()
1k k k k x
x
d
α+=+
产生的迭代序列在无舍入误差的假定下，最多经过n 次迭代就能求得()f x 的最小点，也就是方程组A x
b =的解。
共轭梯度法中关键的两点是，确定迭代格式中的搜索方向和最佳步长。搜索方向()k
d ，与前一次的搜索方向关于(
)
1k d -关于矩阵A 共轭，即
(
)
()
()1,0k k d A d
+=，
然后从点
()
k x
出发，沿
()
k d
方向求得
()f
x 的极小值点
()
1k x
+，即
(
)
()
()(
)
()
()0
m in k k k
k f
x d
f x d
ααα≥+=+ 。
由此解得最佳步长k α和参数k β的表达式为
()()()()
()()()()
1, k T
k k T k k
k
k T
k k T
k r
d
r
A d d A d
d
A d
α
β
+==-
共轭梯度法的计算公式为：
()()()
()()()()
()()()()()()
()()
()()00011(1)
11(1)
r r
r k k T k k k T k k k k k k k T k k k T k k k k k d b A x r d d A d
x x d b A x r
A d d A d d
d ααββ++++++?==-?
?
=??
?=+??=-???=-??=+??
程序框图
计算实习
用共轭梯度法求解线性方程组x
A b
=，其中
21112
1
, 12101
21A b --????????
-?????
???==????-????
????--?
???
矩阵A 的阶数n 分别取为100,200,400，指出计算结果是否可靠。 Matlab 代码
%使用说明：共轭梯度法求解Ax=b %命令行中输入矩阵A 及b %然后调用函数getd 进行计算
%变量含义：n —方程阶数，x0—初始向量 %e —计算精度，r —残向量 clear clc format long
n=input('请输入方程阶数n='); %输入矩阵阶数 A=zeros(n); b=zeros(n,1);
A(1,1)=-2;A(1,2)=1;A(n,n-1)=1;A(n,n)=-2; b(1,1)=-1;b(n,1)=-1; for i=2:n-1;
A(i,i)=-2;A(i,i-1)=1;A(i,i+1)=1; end; A;
b; %生成对应阶数的矩阵A 和b x0=zeros(n,1); %生成初始向量 e=input('请输入计算精度e='); %输入计算精度 x=getd(A,b,x0,e); %调用函数getd 进行计算
if n==100 %对x元素进行重新排列
x=reshape(x,10,10)
elseif n==200
x=reshape(x,10,20)
else
x=reshape(x,20,20)
end
被调用函数getd
function x=getd(A,b,x0,e)% 矩阵A,b，初始向量x0，精度e n=size(A,1);% 获取矩阵A的阶数
x=x0;%初始向量
r=b-A*x;%残向量
d=r;%搜索向量
for k=0:(n-1)
p=(r'*r)/(d'*A*d);
x=x+p*d;
r2=b-A*x;
if ((norm(r2)<=e)||(k==n-1))
x;
break;
end
q=norm(r2)^2/norm(r)^2;
d=r2+q*d;
r=r2;
end
4.三次样条插值计算原理
程序框图
计算实习 给定函数
2
1(x )(1x 1)
115f x
=
-≤≤+ .取等距节点，构造三次样条插值10(x )S 。
Matlab 代码
%使用说明：该程序解决的是三次样条插值中，第1,2类边界条件的问题 %各变量含义：a,b —插值区间左右端点 % n —插值节点数目 % p,q —左右端点导数值
% A ，M ，d —用于求解AM=d 中，矩阵M 的值 % C —存放各区间内插值函数的系数矩阵
% zglu —利用追赶法进行LU 分解，求解AM=d 的函数 clear clc format long
%输入区间，计算插值节点 a=input('请输入区间左端点a='); b=input('请输入区间右端点b=');
n=input('请输入插值节点数目(包括左右端点)n='); d=zeros(n,1);x=zeros(n,1);y=zeros(n,1);A=zeros(n); h=(b-a)/(n-1); fprintf('步长h=%d',h) for i=1:n x(i,:)= a+h*(i-1);
y(i,:)=1/(1+25*(x(i,:))^2);%计算插值节点处的函数值 end
%选择边界条件进行计算,并输入区间左右端点的导数值p ，q xz=input('请选择边界条件类型（1或2或3）xz='); fprintf('以第%d 类边界条件进行计算',xz);
p=input('请输入区间左端边界条件p=');
q=input('请输入区间右端边界条件q=');
%计算矩阵A及矩阵d
if xz==1
A(1,1)=2;A(1,2)=1/2;A(n,n)=2;A(n,n-1)=1/2;
for j=1:n-1
d(j,:)=(3/h)*((y(j+1,:)-y(j,:))/h-(y(j,:)-y(j-1,:))/h);
A(j,j)=2;A(j,j-1)=0.5;A(j,j+1)=0.5;
end
d(1,:)=d(1,:)-1/2*p;
d(n,:)=d(n,:)-1/2*q
elseif xz==2
d(1,:)=(6/h)*((y(2,:)-y(1,:))/h-p);
d(n,:)=(6/h)*(q-(y(n,:)-y(n-1,:))/h);
A(1,1)=2;A(1,2)=1;A(n,n)=2;A(n,n-1)=1;
for j=2:n-1
d(j,:)=(3/h)*((y(j+1,:)-y(j,:))/h-(y(j,:)-y(j-1,:))/h);
A(j,j)=2;A(j,j-1)=0.5;A(j,j+1)=0.5;
end
end
%调用函数zglu用追赶法计算AM=d
M=zglu(A,d);
%各插值区间内函数表达式，系数矩阵为n*4阶矩阵C
for k=2:n
C(k-1,1)=(M(k,:)-M(k-1,:))/(6*h);
C(k-1,2)=(x(k,:)*M(k-1,:)-x(k-1,:)*M(k,:))/(2*h);
C(k-1,3)=1/(2*h)*(x(k-1,:)^2*M(k,:)-x(k,:)^2*M(k-1,:))+1/h*(y(k,:)-1/6*h^2*M(k,:)-
y(k-1,:)-1/6*h^2*M(k-1,:));
C(k-1,4)=1/(6*h)*(x(k,:)^3*M(k-1,:)-x(k-1,:)^3*M(k,:))+1/h*(x(k,:)*(y(k-1,:)-h^2/6* M(k-1,:))-x(k-1,:)*(y(k,:)-h^2/6*M(k,:)));
end
%显示输入数据
disp('您输入的数据如下：')
disp('插值节点x：')
x(:,:)
disp('插值节点y：')
y(:,:)
disp('计算得到矩阵M：')
M(:,:)
%输出插值函数S(x)的表达式
disp('S(x)的表达式为：')
for l=1:n-1
disp([num2str(C(l,1)),'x^3+',num2str(C(l,2)),'x^2+',num2str(C(l,3)),'x+',num2str(C(l,4 )),' ',num2str(x(l,:)),'≤x≤',num2str(x(l+1,:))]);
end
被调用函数zglu
% 追赶法求解三对角方程组
function x=zglu(A,b)
[~,n]=size(A);
L=eye(n);U=zeros(n);y=zeros(n,1);x=zeros(n,1);
U(1,1)=A(1,1);y(1,1)=b(1,1);
计算方法上机实验报告
. / 《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：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; end x=double(x1) K 四、运行结果：
数值分析上机作业
数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：
课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )
则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
东南大学数值分析上机题答案
数值分析上机题 第一章 17.（上机题）舍入误差与有效数 设∑=-= N j N j S 2 2 11 ，其精确值为）111-23（21+-N N 。 （1）编制按从大到小的顺序1 -1 ···1-311-21222N S N +++=，计算N S 的通用 程序； （2）编制按从小到大的顺序1 21 ···1)1(111 222-++--+ -=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时用单精度）； （4）通过本上机题，你明白了什么？ 解： 程序： （1）从大到小的顺序计算1 -1 ···1-311-21222N S N +++= ： function sn1=fromlarge(n) %从大到小计算sn1 format long ; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end （2）从小到大计算1 21 ···1)1(111 2 22 -++--+-= N N S N function sn2=fromsmall(n) %从小到大计算sn2 format long ; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end （3） 总的编程程序为： function p203()
clear all format long; n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%f\n',sn); sn1=fromlarge(n); fprintf('从大到小计算的值为%f\n',sn1); sn2=fromsmall(n); fprintf('从小到大计算的值为%f\n',sn2); function sn1=fromlarge(n) %从大到小计算sn1 format long; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end function sn2=fromsmall(n) %从小到大计算sn2 format long; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end end 运行结果：
计算方法上机题答案
2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量 （1）二分法 （局部，大图不太看得清，故后面两小题都用局部截图） （2）迭代法
（3）牛顿法 顺序消元法 #include for(k=p+1;k 计算方法上机报告 姓名： 学号： 班级：能动上课班级： 题目及求解： 一、对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=0681581482184161n n n n S n ，要求： ① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算； 1 算法思想 （1）根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114 16818485861681 n n n a n n n n n ε??= ---<< ?+++++??； （2）为了保证计算结果的准确性，写程序时，从后向前计算； （3）使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数) 2 算法结构 ;0=s ?? ? ??+-+-+-+= 681581482184161n n n n t n ； for 0,1,2,,n i =??? if 10m t -≤ end; for ,1,2,,0n i i i =--??? ;s s t =+ 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=7，此时求解得： s =3.1415926536； 若保留30位有效数字时，则n=22, 此时求解得： s =3.8。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 二、某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示： 第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1，其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________') 三、结果 从结果可以看出有效位数是6位。 感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。 数值分析上机题目 1、 分别用不动点迭代与Newton 法求解方程250x x e -+=的正根与负根。 2、 Use each of the following methods to find a solution in [0.1,1] accurate to within 10^-4 for 4326005502002010x x x x -+--= a. Bisection method b. Newton’s method c. Secant method d. Method of False Position e. Muller’s method 3、 应用Newton 法求f （x ）的零点，e=10^-6,这里f （x ）=x-sin （x ）。 再用求重根的两种方法求f （x ）的零点。 4、 应用Newton 法求f （x ）的零点，e=10^-6,f(x)=x-sin(x) 再用Steffensen’s method 加速其收敛。 5、 用Neville’s 迭代差值算法，对于函数2 1 (),11125f x x x = -≤≤+进行lagrange 插值。取不同的等分数n=5,10，将区间[-1,1]n 等分，取等距节点。把f(x)和插值多项式的曲线画在同一张图上进行比较。 6、 画狗的轮廓图 7、 Use Romberg integration to compute the following approximations to ? a 、 Determine R1,1,R2,1,R3,1,R4,1and R5,1,and use these approximations to predict the value of the integral. b 、 Determine R2,2 ,R3,3 ,R4,4 ,and R5,5,and modify your prediction. c 、 Determine R6,1 ,R6,2 ,R6,3 ,R6,4 ,R6,5 and R6,6,and modify your prediction. 数值计算方法上机实习题 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 ，并用牛顿迭代法； 《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：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; end x=double(x1) K 四、运行结果： 实验二： 数值计算方法I 上机实验考试题（两题任选一题） １．小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. A. 建立火箭升空过程的数学模型（微分方程）； B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． ２．小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为 0)0(,0,01222==≤≤-+?? ? ??-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻． 今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））： 现有两种估计比例系数k 的方法： 1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。 2．用这组数据拟合一个k ． 请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）． 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0, 解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志 ,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n 数值分析上机题1 设2 21 1N N j S j ==-∑ ，其精确值为1311221N N ??-- ?+?? 。 （1）编制按从大到小的顺序222 111 21311 N S N = +++---，计算N S 的通用程序。 （2）编制按从小到大的顺序22 21111(1)121 N S N N =+++----，计算N S 的通用程序。 （3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数。（编制程序时用单精度） （4）通过本上机题，你明白了什么？ 程序代码（matlab 编程）： clc clear a=single(1./([2:10^7].^2-1)); S1(1)=single(0); S1(2)=1/(2^2-1); for N=3:10^2 S1(N)=a(1); for i=2:N-1 S1(N)=S1(N)+a(i); end end S2(1)=single(0); S2(2)=1/(2^2-1); for N=3:10^2 S2(N)=a(N-1); for i=linspace(N-2,1,N-2) S2(N)=S2(N)+a(i); end end S1表示按从大到小的顺序的S N S2表示按从小到大的顺序的S N 计算结果 通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。从大到小的顺序计算得到的结果的有效位数少。计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。 实习题1 1用两种不容的顺序计算644834.11000 12≈∑=-n n ，分析误差的变化 （1）顺序计算 源代码： #include if(n%1000==0) printf("sum[%d]=%-30f",n,sum); if(n<=1)break; n--; } printf("sum[%d]=%f\n",n,sum); } 结果： 2已知连分数 )) / /(... /( 3 2 2 1 1 n n b a a b a b a b f + + + + = 利用下面的方法计算f: 1 1)0 ,..., 2 ,1 ( , d f n n i d a b d b d i i i i n n = - - = + = = + + 写一个程序，读入n, n n b a,,计算并打印f 源代码： #include 本科实验报告 课程名称：计算机数值方法 实验项目：方程求根、线性方程组的直接解 法、线性方程组的迭代解法、代数插值和最 小二乘拟合多项式 实验地点：行勉楼 专业班级： ******** 学号： ********* 学生姓名： ******** 指导教师：李誌，崔冬华 2016年 4 月 8 日 y = x*x*x + 4 * x*x - 10; return y; } float Calculate(float a,float b) { c = (a + b) / 2; n++; if (GetY(c) == 0 || ((b - a) / 2) < 0.000005) { cout << c <<"为方程的解"<< endl; return 0; } if (GetY(a)*GetY(c) < 0) { return Calculate(a,c); } if (GetY(c)*GetY(b)< 0) { return Calculate(c,b); } } }; int main() { cout << "方程组为：f(x)=x^3+4x^2-10=0" << endl; float a, b; Text text; text.Getab(); a = text.a; b = text.b; text.Calculate(a, b); return 0; } 2.割线法： // 方程求根（割线法）.cpp : 定义控制台应用程序的入口点。// #include "stdafx.h" #include"iostream" 心得体会 使用不同的方法，可以不同程度的求得方程的解，通过二分法计算的程序实现更加了解二分法的特点，二分法过程简单，程序容易实现，但该方法收敛比较慢一般用于求根的初始近似值，不同的方法速度不同。面对一个复杂的问题，要学会简化处理步骤，分步骤一点一点的循序处理，只有这样，才能高效的解决一个复杂问题。 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。 Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear %%%给定求解区间 b=1.5; a=0; %%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f); 计算方法实验报告格式 小组名称: 组长姓名(班号)： 小组成员姓名(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一个完整的实验，应包括数据准备、理论基础、实验内容及方法，最终对实验结果进行分析，以达到对理论知识的感性认识，进一步加深对相关算法的理解，数值实验以实验报告形式完成，实验报告格式如下： 一、实验名称 实验者可根据报告形式需要适当写出． 二、实验目的及要求 首先要求做实验者明确，为什么要做某个实验，实验目的是什么，做完该实验应达到什么结果，在实验过程中的注意事项，实验方法对结果的影响也可以以实验目的的形式列出． 三、算法描述（实验原理与基础理论） 数值实验本身就是为了加深对基础理论及方法的理解而设置的，所以要求将实验涉及到的理论基础，算法原理详尽列出． 四、实验内容 实验内容主要包括实验的实施方案、步骤、实验数据准备、实验的算法以及可能用到的仪器设备． 五、程序流程图 画出程序实现过程的流程图，以便更好的对程序执行的过程有清楚的认识，在程序调试过程中更容易发现问题． 六、实验结果 实验结果应包括实验的原始数据、中间结果及实验的最终结果，复杂的结果可以用表格 形式列出，较为简单的结果可以与实验结果分析合并出现． 七、实验结果分析 实验结果分析包括对对算法的理解与分析、改进与建议． 数值实验报告范例 为了更好地做好数值实验并写出规范的数值实验报告，下面给出一简单范例供读者参考． 数值实验报告 小组名称: 小组成员(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一、实验名称 误差传播与算法稳定性． 二、实验目的 1．理解数值计算稳定性的概念． 2．了解数值计算方法的必要性． 3．体会数值计算的收敛性与收敛速度． 三、实验内容 计算dx x x I n n ? += 1 10 ，1,2,,10n = ． 四、算法描述 由 dx x x I n n ? += 1 10 ，知 dx x x I n n ?+=--101110，则 如有帮助欢迎下载支持 数值分析上机题 姓名：陈作添 学号： 040816 习题 1 20．（上机题）舍入误差与有效数 N 1 1 3 1 1 设 S N ，其精确值为 。 2 2 2 N N 1 j 2 j 1 （1）编制按从大到小的顺序 1 1 1 ，计算 S 的通用程序。 S N 1 32 1 N 2 1 N 2 2 （2）编制按从小到大的顺序 1 1 1 ，计算 S 的通用程序。 S N 1 (N 1)2 1 22 1 N N 2 （3）按两种顺序分别计算 S 102 ， S 104 ， S 106 ，并指出有效位数。 （编制程序时用单精度） （4）通过本上机题，你明白了什么？ 按从大到小的顺序计算 S N 的通用程序为： 按从小到大的顺序计算 S N 的通用程序为： #include 上海电力学院 数值分析上机实验报告 题目：数值分析上机实验报告 学生姓名：11111111111 学号：111111********* 专业：1111 2013年12月30日 数值计算方法上机实习题 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+ -=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用n I I n n 51 511+-=-，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 (1) 解答：n=0， 0.1823)05ln()15ln()5(5151510101 0=+-+=++=+=+=???x d x dx x dx x x I n n 这里可以用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 >> I I20=-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<21 End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> I I =0.0083 (3) 算法误差分析： 计算在递推过程中传递截断误差和舍入误差 第一种算法：（从1——>20） * 000 e I I =- * **21111120 11 5(5)5()555n n n n n n n n n n e I I I I I I e e e n n ------=-=-+ --+=-=== 1 x x 205x +????? d 7.99810 3 -?= 计算方法B 上机报告 第1题 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示： (1)请用合适的曲线拟合所测数据点； (2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图； 问题分析和算法思想： 本题的主要目的是对21个测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值，可以用的插值方法很多：多项式插值、Lagrange 插值、Newton 插值、三次样条插值等。由于数值点较多时，采用高次多项式插值将产生很大的误差，用拉格朗日插值多项式会出现龙格现象。故为了将所有的数据点都用上，且题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值的方法较为合适。 计算光缆长度近似值，只需将每两点之间的距离算出，然后依次相加，所得的折线长度，即为光缆长度的近似值。 光缆长度计算公式： 19 1 k k k l +===∑? ? ? 算法结构： 三次样条算法结构见《计算方法教程》P110。 源程序： clear;clc; x=0: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]; d=y; plot(x,y,'k.','markersize',15) hold on %%%计算二阶差商 for k=1:2 for i=21:-1:(k+1) d(i)=(d(i)-d(i-1))/(x(i)-x(i-k)); end end %%%假定d的边界条件，采用自然三次样条 for i=2:20 d(i)=6*d(i+1); end d(1)=0; d(21)=0; %%%追赶法求解带状矩阵的m值 a=0.5*ones(1,21); b=2*ones(1,21); c=0.5*ones(1,21); a(1)=0;c(21)=0; u=ones(1,21); u(1)=b(1); r=c; yy(1)=d(1); %%%追的过程 for k=2:21 l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); yy(k)=d(k)-l(k)*yy(k-1); end %%%赶的过程 m(21)=yy(21)/u(21); for k=20:-1:1 m(k)=(yy(k)-r(k)*m(k+1))/u(k); end %%%利用插值点画出拟合曲线 k=1; nn=100; xx=linspace(0,20,nn); l=0; for j=1:nn for i=2:20 if xx(j)<=x(i) k=i;西安交通大学计算方法B上机报告
数值分析上机题目详解
数值分析上机题目
数值计算方法上机实习题
计算方法上机实验报告
数值计算方法I上机实验考试题
(完整版)数值计算方法上机实习题答案
《数值计算方法》上机实验报告
东南大学《数值分析》-上机题
计算方法与实习上机题答案
太原理工大学数值计算方法实验报告
(完整版)哈工大-数值分析上机实验报告
计算方法实验报告格式
数值分析上机题参考答案.docx
计算方法上机答案
计算方法B上机报告