东南大学数值分析上机作业汇总
数值分析大作业
数值分析上机作业（一）一、算法的设计方案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 。
数值分析报告上机题课后作业全部-东南大学
实用标准文案文档大全上机作业题报告2015.1.9 USER1.Chapter 11.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） （4）通过本次上机题，你明白了什么？1.2程序1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。
按从小到大的顺序，有效位数分别为：5，6，6。
可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。
因此，采取从小到大的顺序累加得到的结果更加精确。
2.Chapter 22.1题目（1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。
（2）给定方程03)(3=-=x xx f ,易知其有三个根3,0,3321=*=*-=*x x x○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。
试确定尽可能大的δ。
○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。
算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。
运行Find.m，得到在不同的搜索精度下的最大sigma值。
东南大学数值分析上机题答案
东南⼤学数值分析上机题答案数值分析上机题第⼀章17.（上机题）舍⼊误差与有效数设∑=-=Nj N j S 2211，其精确值为）111-23（21+-N N 。
（1）编制按从⼤到⼩的顺序1-1···1-311-21222N S N +++=，计算N S 的通⽤程序；（2）编制按从⼩到⼤的顺序121···1)1(111222-++--+-=N N S N ，计算NS 的通⽤程序；（3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时⽤单精度）；（4）通过本上机题，你明⽩了什么？解：程序：（1）从⼤到⼩的顺序计算1-1···1-311-21222N S N +++=：function sn1=fromlarge(n) %从⼤到⼩计算sn1format long ; sn1=single(0); for m=2:1:nsn1=sn1+1/(m^2-1); end end（2）从⼩到⼤计算121···1)1(111222-++--+-=N N S N function sn2=fromsmall(n) %从⼩到⼤计算sn2format long ; sn2=single(0); for m=n:-1:2sn2=sn2+1/(m^2-1); end end（3）总的编程程序为： function p203()clear allformat 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:nsn1=sn1+1/(m^2-1);endendfunction sn2=fromsmall(n) %从⼩到⼤计算sn2 format long;sn2=single(0);for m=n:-1:2sn2=sn2+1/(m^2-1);endendend运⾏结果：从⽽可以得到N值真值顺序值有效位数2 100.740050 从⼤到⼩0.740049 5从⼩到⼤0.740050 64 100.749900 从⼤到⼩0.749852 3从⼩到⼤0.749900 66 100.749999 从⼤到⼩0.749852 3从⼩到⼤0.749999 6（4）感想：通过本上机题，我明⽩了，从⼩到⼤计算数值的精确位数⽐较⾼⽽且与真值较为接近，⽽从⼤到⼩计算数值的精确位数⽐较低。
东南大学数值分析上机
第一章一、题目设∑=-=Nj N j S 2211，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算SN 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算SN 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 二、MATLAB 程序N=input('请输入N(N>1)：');AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); %single 使其为单精度Sn1=single(0); %从小到大的顺序 for a=2:N; Sn1=Sn1+1/(a^2-1); endSn2=single(0); %从大到小的顺序 for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); endfprintf('Sn 的值 (N=%d)\n',N);disp('____________________________________________________') fprintf('精确值 %f\n',AccurateValue); fprintf('从大到小计算的结果 %f\n',Sn1);fprintf('从小到大计算的结果 %f\n',Sn2);disp('____________________________________________________')三、结果请输入N(N>1)：100Sn的值 (N=100)____________________________________________________精确值 0.740049从大到小计算的结果 0.740049从小到大计算的结果 0.740050____________________________________________________请输入N(N>1)：10000Sn的值 (N=10000)____________________________________________________精确值 0.749900从大到小计算的结果 0.749852从小到大计算的结果 0.749900____________________________________________________请输入N(N>1)：1000000Sn的值 (N=1000000)____________________________________________________精确值 0.749999从大到小计算的结果 0.749852从小到大计算的结果 0.749999____________________________________________________四、结果分析可以得出，算法对误差的传播又一定的影响，在计算时选一种好的算法可以使结果更为精确。
东南大学出版社第二版《数值分析》上机作业答案(前三章)
for (i=k+1;i<N;i++) // { lik=a[i][k]/a[k][k]; //实施消去过程，得到上三角系数增广矩阵 for (j=k;j<M;j++) // { a[i][j]=a[i][j]‐lik*a[k][j]; // } } } cout<<"经列主元高斯消去法得到的数组为:"<<endl; // for (b=0;b<N;b++) // { cout<<endl; //输出经过列主元消去法处理过的系数增广矩阵 for (c=0;c<M;c++) { cout<<setw(7)<<a[b][c]; // } } cout<<endl; double x[N]; // double s; int f,g; x[N‐1]=a[N‐1][M‐1]/a[N‐1][N‐1]; // for (f=N‐2;f>=0;f‐‐) // { s=0; for (g=f+1;g<N;g++) //由上三角形的系数增广矩阵求出方程组的解 { s=s+a[f][g]*x[g]; // } x[f]=(a[f][N]‐s)/a[f][f]; // } cout<<"方程组的解为:"<<endl; for (b=0;b<N;b++) //输出方程组的解 {
1
当 n=10000 时，s3=0.7499 Press any key to continue （分析 S1 的 6 位数字中，有效位数为 4 位； S2 的所有数字都是有效数字。 ） 当 n=1000000 时，s1=‐14.2546 当 n=1000000 时，s2=‐14.2551 当 n=1000000 时，s3=0.749999 Press any key to continue （分析： S1 的 6 位数字中，没有有效数字； S2 的 6 位数字中，没有有效数字。 ） 由运行结果可知，当精度比较低时，按从大数开始累加到小数的计算结果的精度低于按从小数 累加到大数的计算结果的精度。 至于当 n=1000000 时，S1 和 S2 得出了负数结果，可能是由于循环次数过多，导致数据溢出， 从而得出错误结果。 习题 2 20．程序如下： //给定误差限为：0.5e‐6 //经过试算得当 delta 最大取道 0.7745966 时，迭代得到的根都收敛于 0 #include <iostream.h> #include <math.h> void main () { double x,u; int count=0; u=10.0; cout<<"请输入 x 的初值"<<endl; cin>>x; for (count=0;abs(u)>5;count++) { x=x‐(x*x*x‐3*x)/(3*(x*x‐1)); u=10000000*x; if(count>5000) { cout<<"迭代结果不收敛于 0！"<<endl; break; } } cout<<"x="<<x<<endl<<endl;
数值分析上机题Matlab(东南大学)3
0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72
152 139 128 119 110 103 96 90 85 80 76 72 68 65 62 59 56 53 51 49 47 45 43 41 39 38
========================================================================================================================
======================================================================================================================================================================== 习题 3_36 ======================================================================================================================================================================== Omega n x1 x2 x3 x4 x5 x6 x7 x8 x9
-0.71279 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281
东南大学数值分析--精选上机作业汇总--精选.doc
在x0∈（-1，
）区间，取以下初值，分别调用
newton.m函数文件，得到
结果如下：
X0
X1
迭代次数
-0.95
1.732051
9
-0.85
1.732051
6
-0.80
1.732051
10
-0.78
1.732051
0.0000005
由newton1.m的运行过程表明，在整个区间上均收敛于0，即根x*2。
在x0∈（，1）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：
X0
X1
迭代次数
0.80
-1.732051
10
0.90
-1.732051
7
0.95
-1.732051
9
0.98
-1.732051
15
计算结果显示，迭代序列局部收敛于1.730251，即根x3*。
在x0∈（
，
）区间，取以下初值，分别调用
newton.m函数文件，得到
结果如下：
X0
X1
迭代次数
-0.70
0.000000
5
-0பைடு நூலகம்20
0.000000
3
-0.05
0.000000
3
4
0.05
0.20
0.70
0.0000003
0.0000003
1
32
1
N2
1
（2）编制按从大到小的顺序
SN
1
1
1
，计算SN的通用
数值分析上机作业
数值分析上机作业2实验一：（1） ①用不动点迭代法求()013=--=x x x f 的根，至少设计两种迭代格式，一个收敛一个发散，1210-=ε。
（2） ②对迭代格式使用Aitken 加速，观察敛散性变化。
1取递推公式31)1(+=x x ，可以得到收敛时的迭代结果为：x=(2)^(1/3); t=1; while(1)if(abs(x-(x+1)^(1/3))<10^-12) break; endx=(x+1)^(1/3);t=t+1;end t xtemp=x^3-x-1 %带回来验证下 t = 16 x =1.324717957243755 temp =-4.225952920933196e-12 加速后代码如下 x=1;x=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)) t=1; while(1)if(abs(x-(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)))<10^-10) break; endx=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)); t=t+1; if (t>100000)break; %防止运算次数过多 end endfprintf('需要%d 次',t);输出需要670次>>此处取10^10-=ε，若到10^-12次方则可能需要运行更多次取13-=x x 则迭代发散。
使用aitken 加速计算结果如下 x=2; t=1; while(1)if( abs(x-((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1))<10^-10) break; end t=t+1;x=((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1); if (t>100000)break; %防止运算次数过多 end end t xt = 108x = 1.324717956244172由此可见经过aitken 加速以后，原来发散的迭代格式收敛了。
(完整word版)数值分析上机作业1-1解析
考虑一个高次的代数多项式
（E1-1）
显然该多项式的全部根为l，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。现考虑该多项式方程的一个扰动
(E1-2)
其中 是一个非常小的数。这相当于是对（E1-1）中 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。
ve=zeros(1,21);
ve(21-Numb)=ess;
root=roots(poly(1:20)+ve);
x0=real(root); y0=imag(root);
plot(x0',y0','*');
disp(['对扰动项 ',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为:']);
ess分别为1e-6,1e-8.1e-10,1e-12的图像如下：
从实验的图形中可以看出，当ess充分小时，方程E.1.1和方程E.1.2的解相差很小，当ess逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。
（2）将扰动项加到x18上后，ess=1e-009时方程的解都比较准确，没有出现复共轭现象。ess=1e-008时误差与x19(ess=1e-009)时相当，即扰动加到x18上比加到x19小一个数量级。对x8的扰动ess=1000时没有出现复共轭，误差很小；对x的扰动ess=10e10时没有出现复共轭，误差很小。因此，扰动作用到xn上时，n越小，扰动引起的误差越小。
if((Numb>20)|(Numb<0))errordlg('请输入正确的扰动项:[0 20]之间的整数!');return;end
东南大学数值分析上机
第一章一、题目设∑=-=Nj N j S 2211，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算SN 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算SN 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 二、MATLAB 程序N=input('请输入N(N>1)：');AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); %single 使其为单精度 Sn1=single(0); %从小到大的顺序 for a=2:N; Sn1=Sn1+1/(a^2-1); endSn2=single(0); %从大到小的顺序 for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); endfprintf('Sn 的值 (N=%d)\n',N);disp('____________________________________________________') fprintf('精确值 %f\n',AccurateValue); fprintf('从大到小计算的结果 %f\n',Sn1); fprintf('从小到大计算的结果 %f\n',Sn2);disp('____________________________________________________')三、结果请输入N(N>1)：100Sn的值(N=100)____________________________________________________精确值0.740049从大到小计算的结果0.740049从小到大计算的结果0.740050____________________________________________________请输入N(N>1)：10000Sn的值(N=10000)____________________________________________________精确值0.749900从大到小计算的结果0.749852从小到大计算的结果0.749900____________________________________________________请输入N(N>1)：1000000Sn的值(N=1000000)____________________________________________________精确值0.749999从大到小计算的结果0.749852从小到大计算的结果0.749999____________________________________________________四、结果分析可以得出，算法对误差的传播又一定的影响，在计算时选一种好的算法可以使结果更为精确。
东南大学数值分析上机报告完整版
数值分析上机实验报告目录1.chapter1舍入误差及有效数 (2)2.chapter2Newton迭代法 (3)3.chapter3线性代数方程组数值解法-列主元Gauss消去法 (7)4.chapter3线性代数方程组数值解法-逐次超松弛迭代法 (9)5.chapter4多项式插值与函数最佳逼近 (10)1.chapter1舍入误差及有效数1.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度）（4）通过本次上机题，你明白了什么？1.2编写相应的matlab 程序clear;N=input('please input N:');AValue=((3/2-1/N-1/(N+1))/2);sn1=single(0);sn2=single(0);for i=2:Nsn1=sn1+1/(i*i-1); %从大到小相加的通用程序%endep1=abs(sn1-AValue);for j=N:-1:2sn2=sn2+1/(j*j-1); %从小到大相加的通用程序%endep2=abs(sn2-AValue);fprintf('精确值为:%f\n',AValue);fprintf('从大到小的顺序累加得sn=%f\n',sn1);fprintf('从大到小相加的误差ep1=%f\n',ep1);fprintf('从小到大的顺序累加得sn=%f\n',sn2);fprintf('从小到大相加的误差ep2=%f\n',ep2);disp('=================================');1.3matlab 运行程序结果>> chaper1please input N:100精确值为:0.740050从大到小的顺序累加得sn=0.740049从大到小相加的误差ep1=0.000001从小到大的顺序累加得sn=0.740050从小到大相加的误差ep2=0.000000>> chaper1please input N:10000精确值为:0.749900从大到小的顺序累加得sn=0.749852从大到小相加的误差ep1=0.000048从小到大的顺序累加得sn=0.749900从小到大相加的误差ep2=0.000000>> chaper1please input N:1000000精确值为:0.749999从大到小的顺序累加得sn=0.749852从大到小相加的误差ep1=0.000147从小到大的顺序累加得sn=0.749999从小到大相加的误差ep2=0.0000001.4结果分析以及感悟按照从大到小顺序相加的有效位数为：5,4,3。
数值分析上机作业(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,用试算法确定最佳松弛因子。
东南大学数值分析上机题上.doc
数值分析上机报告姓名：学号：专业：2013年10月27日第一章舍入误差与有效数 设2211NN j S j ==-∑，其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭。
（1）编制按从大到小的顺序22211121311N S N =+++---，计算N S 的通用程序。
（2）编制按从小到大的顺序2221111(1)121N S N N =+++----，计算NS 的通用程序。
（3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数。
（编制程序时用单精度） （4）通过本上机题，你明白了什么？解： （1）#include<stdio.h> void main() { float n,i,s; printf("please input n="); scanf("%f",&n); for(i=2,s=0;i<=n;s+=1/(i*i-1),i++); printf("s=%f\n",s); }（2）#include<stdio.h> void main() { float n,i,s; printf("please input n="); scanf("%f",&n); for(i=n,s=0;i>=2;s+=1/(i*i-1),i--); printf("s=%f\n",s); }（3）按从大到小顺序：210S =0.740049 有效位数6位 410S =0.749852 有效位数3位 610S =0.749852 有效位数3位 按从小到大顺序：210S =0.740050 有效位数5位410S =0.749900 有效位数6位610S =0.749999 有效位数6位(4)通过上述实验数据可以看出此次算法使用从小到大的顺序进行得到的数据相对而言更精确，可以得到这样的启示：在计算数值时，要先分析不同算法对结果的影响，避免大数吃小数的现象，找出能得到更精确的结果的算法。
东南大学数值分析上机实验题(下)
数值分析上机报告XX：学号：2013年12月22日第四章38.（上机题）3次样条插值函数（1）编制求第一型3次样条插值函数的通用程序;端点条件为'0y =0.8，'10y =0.2。
用所编制程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,…9)。
解：（1）#include <iostream.h> #include <math.h>floatx1[100],f1[100],f2[99],f3[98],m[100],a[100][101],x,d[100]; float c[99],e[99],h[99],u[99],w[99],y_0,y_n,arr ,s; int i,j,k,n,q;void selectprint(float y) {if ((y>0)&&(y!=1)) cout<<"+"<<y; else if (y==1) cout<<"+"; else if (y<0) cout<<y; }void printY(float y){ if (y!=0) cout<<y; }float calculation(float x){ for (j=1;j<=n;j++) if (x<=x1[j]) {s=(float)(f1[j-1]+c[j-1]*(x-x1[j-1])+m[j-1]/2.0*(x-x1[j-1])*(x-x1[j-1])+e[j-1]*(x-x1[j-1])*(x-x1[j-1])*(x-x1[j-1])); break; }return s; }void main() {do{cout<<"请输入n值:";cin>>n;if ((n>100)||(n<1)) cout<<"请重新输入整数(1..100):"<<endl;}while ((n>100)||(n<1));cout<<"请输入Xi(i=0,1,...,"<<n<<"):";for (i=0;i<=n;i++) cin>>x1[i];cout<<endl;cout<<"请输入Yi(i=0,1,...,"<<n<<"n):";for (i=0;i<=n;i++) cin>>f1[i];cout<<endl;cout<<"请输入第一型边界条件Y'0，Y'n:";cin>>y_0>>y_n;cout<<endl;for (i=0;i<n;i++) h[i]=x1[i+1]-x1[i];for (i=1;i<n;i++) u[i]=(float) (h[i-1]/(h[i-1]+h[i]));for (i=1;i<n;i++) w[i]=(float) (1.0-u[i]);for (i=0;i<n;i++) f2[i]=(f1[i+1]-f1[i])/h[i]; //一阶差商for (i=0;i<n-1;i++) f3[i]=(f2[i+1]-f2[i])/(x1[i+2]-x1[i]); //二阶差商for (i=1;i<n;i++) d[i]=6*f3[i-1]; //求出所有的d[i]d[0]=6*(f2[0]-y_0)/h[0];d[n]=6*(y_n-f2[n-1])/h[n-1];for (i=0;i<=n;i++)for (j=0;j<=n;j++)if (i==j) a[i][j]=2;else a[i][j]=0;a[0][1]=1;a[n][n-1]=1;for (i=1;i<n;i++){a[i][i-1]=u[i];a[i][i+1]=w[i];}for (i=0;i<=n;i++) a[i][n+1]=d[i];for (i=1;i<=n;i++) //用追赶法解方程,得M[i]{arr=a[i][i-1];for (j=0;j<=n+1;j++)a[i][j]=a[i][j]-a[i-1][j]*arr/a[i-1][i-1];}m[n]=a[n][n+1]/a[n][n];for (i=n-1;i>=0;i--) m[i]=(a[i][n+1]-a[i][i+1]*m[i+1])/a[i][i];for (i=0;i<n;i++) //计算S（x）的表达式c[i]=(float) (f2[i]-(1.0/3.0*m[i]+1.0/6.0*m[i+1])*h[i]);for (i=0;i<n;i++)e[i]=(m[i+1]-m[i])/(6*h[i]);for (i=0;i<n;i++){cout<<"X属于区间["<<x1[i]<<","<<x1[i+1]<<"]时"<<endl<<endl;cout<<"S(x)=";printY(f1[i]);if (c[i]!=0){selectprint(c[i]);cout<<"(x";printY(-x1[i]);cout<<")";}if (m[i]!=0){selectprint((float)(m[i]/2.0));for (q=0;q<2;q++){cout<<"(x";printY(-x1[i]);cout<<")";}}if (e[i]!=0){selectprint(e[i]);for (q=0;q<3;q++){cout<<"(x";printY(-x1[i]);cout<<")";}}cout<<endl<<endl;}cout<<"针对本题，计算S(i+0.5)(i=0,1,...,9):"<<endl;for (i=0;i<10;i++){if ((i+0.5<=x1[n])&&(i+0.5>=x1[0])){calculation((float)(i+0.5));cout<<"S("<<i+0.5<<")="<<s<<endl;}else cout<<i+0.5<<"超出定义域"<<endl;};cout<<endl;}（2）编制的程序求车门的3次样条插值函数S(x)：x属于区间[0,1]时；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)x属于区间[1,2]时；S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1) x属于区间[2,3]时；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2)x属于区间[3,4]时；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3) x属于区间[4,5]时；S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4) x属于区间[5,6]时；S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5) x属于区间[6,7]时；S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6) x属于区间[7,8]时；S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7) x属于区间[8,9]时；S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8) x属于区间[9,10]时；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9)S(0.5)=2.90856 S(1.5)=3.67843 S (2.5)=4.38147 S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976 S(9.5)=5.7323第六章21．（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;（5）对于初值问题h=,应用（1）~（4）中的四种方法进行计算，并将计算结果和精确解取步长0.13=+作比较；y x x()3/(1)（6）通过本上机题，你能得到哪些结论？解：#include<iostream.h>#include<fstream.h>#include<stdlib.h>#include<math.h>ofstream outfile("data.txt");//此处定义函数f(x,y)的表达式//用户可以自己设定所需要求得函数表达式double f1(double x,double y){double f1;f1=(-1)*x*x*y*y;return f1;}//此处定义求函数精确解的函数表达式double f2(double x){double f2;f2=3/(1+x*x*x);return f2;}//此处为精确求函数解的通用程序void accurate(double a,double b,double h){double x[100],accurate[100];x[0]=a;int i=0;outfile<<"输出函数准确值的程序结果:\n";do{x[i]=x[0]+i*h;accurate[i]=f2(x[i]);outfile<<"accurate["<<i<<"]="<<accurate[i]<<'\n';i++;}while(i<(b-a)/h+1);}//此处为经典Runge-Kutta公式的通用程序void RK4(double a,double b,double h,double c) {int i=0;double k1,k2,k3,k4;double x[100],y[100];y[0]=c;x[0]=a;outfile<<"输出经典Runge-Kutta公式的程序结果:\n"; do{x[i]=x[0]+i*h;k1=f1(x[i],y[i]);k2=f1((x[i]+h/2),(y[i]+h*k1/2));k3=f1((x[i]+h/2),(y[i]+h*k2/2));k4=f1((x[i]+h),(y[i]+h*k3));y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;outfile<<"y"<<"["<<i<<"]="<<y[i]<<'\n';i++;}while(i<(b-a)/h+1);}//此处为4阶Adams显式方法的通用程序void AB4(double a,double b,double h,double c) {double x[100],y[100],y1[100];double k1,k2,k3,k4;y[0]=c;x[0]=a;outfile<<"输出4阶Adams显式方法的程序结果:\n";for(int i=0;i<=2;i++){x[i]=x[0]+i*h;k1=f1(x[i],y[i]);k2=f1((x[i]+h/2),(y[i]+h*k1/2));k3=f1((x[i]+h/2),(y[i]+h*k2/2));k4=f1((x[i]+h),(y[i]+h*k3));y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}int j=3;y1[0]=y[0];y1[1]=y[1];y1[2]=y[2];y1[3]=y[3];do{x[j]=x[0]+j*h;y1[j+1]=y1[j]+(55*f1(x[j],y1[j])-59*f1(x[j-1],y1[j-1])+37*f1(x[j-2],y1[j-2])-9*f1(x[j-3],y1[ j-3]))*h/24;outfile<<"y1"<<"["<<j<<"]="<<y1[j]<<'\n';j++;}while(j<(b-a)/h+1);}//主函数void main(void){double a,b,h,c;cout<<"输入上下区间、步长和初始值:\n";cin>>a>>b>>h>>c;accurate(a,b,h);RK4(a,b,h,c);AB4(a,b,h,c);}float si(int u,int v){float sum=0; int q;for(q=0;q<k;q++)sum+=a[u][q]*a[q][v];sum=a[u][v]-sum;return sum;}void exchange(int g){int t; float temp;for(t=0;t<n;t++){temp=a[k][t];a[k][t]=a[g][t];a[g][t]=temp;}}void analyze(){int t;float si(int u,int v);for(t=k;t<n;t++)a[k][t]=si(k,t);for(t=(k+1);t<m;t++)a[t][k]=(float)(si(t,k)/a[k][k]);}void ret(){int t,z;float sum;x[m-1]=(float)a[m-1][m]/a[m-1][m-1];for(t=(m-2);t>-1;t--){sum=0;for(z=(t+1);z<m;z++)sum+=a[t][z]*x[z];x[t]=(float)(a[t][m]-sum)/a[t][t];}}（5）由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值y(x i)和精确值和数值解的误差：由AB4方法得出的结果为：Y1[0]=3 y1[1]=2.997 y1[2]=2.97619 y1[3]=2.92113 y1[4]=2.81839 y1[5]=2.66467 y1[6]=2.4652 y1[7]=2.23308 y1[8]=1.98495 y1[9]=1.73704 y1[10]=1.50219 y1[11]=1.28876 y1[12]=1.10072 y1[13]=0.93871 y1[14]=0.801135y1[15]=0.685335（6）本次上机作业让我知道了在遇到复杂问题中，无法给出解析式的情况下，要学会灵活使用各种微分数值解法，并且能计算出不同方法的精度大小。
东南大学数值分析上机报告2
数值分析上机报告作业4一、题目二、算法描述函数S(x)在每个小区间上都是三次多项式，故S’’(x)在小区间上是一次多项式，根据函数值、一阶差值、二阶差值求出d值，再求出M值，回代求出插值函数，最后节点差值输出结果。
三、源程序%求第一型3次样条插值函数的通用程序clearclc% 输入相关参数n = input('Input n: n =');n = n + 1;xn = zeros(1, n);yn = zeros(1, n);xn(1, :) = input('Input x:');yn(1, :) = input('Input y:');% 输入边界条件dy0 = input('Input the derivative of y(0):');dyn = input('Input the derivative of y(n):');% 求d值d = zeros(n, 1);h = zeros(1, n - 1);f1 = zeros(1, n- 1);f2 = zeros(1, n - 2);for i = 1: n - 1h(i) = xn(i + 1) - xn(i); % 一阶差商f1(i) = (yn(i + 1) - yn(i))/h(i);endfor i = 2 : n – 1 % 一二阶差商f2(i) = (f1(i) - f1(i - 1))/(xn(i + 1) - xn(i - 1));d(i) = 6*f2(i);endd(1) = 6*(f1(1) - dy0)/h(1);d(n) = 6*(dyn - f1(n - 1))/h(n - 1);% 求M值A = zeros(n);miu = zeros(1, n -2);lamda = zeros(1, n - 2);for i = 1: n - 2miu(i) = h(i)/(h(i) + h(i + 1));lamda(i) = 1 - miu(i);endA(1, 2) = 1;A(n, n - 1) = 1;for i = 1: nA(i, i) = 2;endfor i = 2: n - 1A(i, i - 1) = miu(i - 1);A(i, i + 1) = lamda(i - 1);endM = A\d;% 回代求插值函数syms x;for i = 1: n - 1Sx(i) = collect(yn(i) + (f1(i) - (M(i)/3 + M(i + 1)/6)*h(i))*(x - ... xn(i)) + M(i)/2*(x - xn(i))^2 + (M(i + 1) - M(i))/(6*h(i))*(x - ... xn(i))^3);Sx(i) = vpa(Sx(i), 4);endS = zeros(1, n - 1);% 求节点插值for i = 1: n - 1x = xn(i) + 0.5;S(i) = yn(i) + (f1(i) - (M(i)/3 + M(i + 1)/6)*h(i)) * (x - xn(i))... + M(i)/2*(x - xn(i))^2 + (M(i + 1) - M(i))/(6*h(i))*(x - xn(i))^3;end% 输出结果disp('S(x) = ');for i = 1: n - 1fprintf(' %s (%d, %d)\n', char(Sx(i)), xn(i), xn(i + 1))disp(' ------------------------------------------------')enddisp('S(i + 0.5)')disp(' i x(i + 0.5) S(i + 0.5)')for i = i: n - 1fprintf(' %d %.4f %.4f\n', i, xn(i) + 0.5, S(i))end四、心得体会使用变量精度算法(VPA)去计算A中每个元素为d小数位精度，其中d是当前设置的位数，结果使每个元素是符号表达式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
东南大学数值分析上机作业汇总-标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII数值分析上机报告院系：学号：姓名：目录作业1、舍入误差与有效数 (1)1、函数文件cxdd.m (1)2、函数文件cddx.m (1)3、两种方法有效位数对比 (1)4、心得 (2)作业2、Newton迭代法 (2)1、通用程序函数文件 (3)2、局部收敛性 (4)（1）最大δ值文件 (4)（2）验证局部收敛性 (4)3、心得 (6)作业3、列主元素Gauss消去法 (7)1、列主元Gauss消去法的通用程序 (7)2、解题中线性方程组 (7)3、心得 (9)作业4、三次样条插值函数 (10)1、第一型三次样条插值函数通用程序： (10)2、数据输入及计算结果 (12)作业1、舍入误差与有效数设∑=-=Nj N j S 2211，其精确值为⎪⎭⎫ ⎝⎛---1112321N N . （1）编制按从小到大的顺序11131121222-⋅⋅⋅+-+-=N S N ，计算N S 的通用程序；（2）编制按从大到小的顺序()12111111222-⋅⋅⋅+--+-=N N S N ，计算N S 的通用程序；（3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序：1、函数文件cxdd.mfunction S=cxdd(N) S=0; i=2.0;while (i<=N) S=S+1.0/(i*i-1); i=i+1;endscript 运行结果（省略>>）：S=cxdd(80) S=0.7375772、函数文件cddx.mfunction S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); endscript 运行结果（省略>>）： S=cddx(80) S=0.7375773、两种方法有效位数对比精确值函数：function S=jqz(N)S=0.5*(1.5-1.0/N-1.0/(N+1));script运行结果（省略>>）4、心得本题重点体现了数值计算中“大数吃小数”的问题，由于计算机计算的截断特点，从大到小的计算会导致小数的有效数被忽略掉。
从题中可以看出，看出按不同的顺序计算的结果是不相同的，按从小到大的顺序计算的值与精确值吻合，而按从大到小的顺序计算的值与精确值有较大的误差。
计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低。
作业2、Newton迭代法（1）给定初值x0及容许误差ε，编制Newton法解方程f(x)=0根的通用程序。
-，x2※=0，x3※=3。
（2）给定方程f(x)=x3/3-x=0,易知其有三个根x1※=3①由Newton方法的局部收敛性可知存在δ＞0，当x0∈（δ-，δ），Newton 迭代序列收敛于根x2※，试确定尽可能大的δ；②试取若干个初始值，观察当x0∈（-∞，-1），（-1，δ-，δ），-），（δ（δ，1），（1，+∞）时，Newton序列是否收敛以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？1、通用程序函数文件定义f(x)函数function f=fun(x)f=x^3/3-x;end定义f(x)导函数function f=dfun(x)f=x*x-1;end定义求近似解函数function [f,n]=newton(x0,ep)flag=1;n=0;while(flag==1)x1=x0-fun(x0)/dfun(x0);n=n+1;if(abs(x1-x0)<=ep||n>100000)flag=0;endx0=x1;endf=x1;endscript运行结果clear;x0=input('请输入初始值x0：');ep=input('请输入容许误差：');[f,n]=newton(x0,ep);fprintf('方程的一个近似解为:%f\n',x1);2、局部收敛性（1）最大δ值文件flag=1;k=1;x0=0;while flag==1sigma=k*10^-6;x0=sigma;k=k+1;m=0;flag1=1;while flag1==1 && m<=10^3x1=x0-fun(x0)/dfun(x0);if abs(x1-x0)<10^-6flag1=0;endm=m+1;x0=x1;endif (flag1==1||abs(x0)>=10^-6)flag=0;endendfprintf('最大值为： %f\n',sigma);运行结果为：最大值为：0.774597x=0的最大区间为即得最大的δ为0.774597，Newton迭代序列收敛于根*2（-0.774597，0.774597）。
（2）验证局部收敛性在x0∈（-∞，-1）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：x。
结果显示，以上初值迭代序列均收敛于-1.732051，即根*1显然，迭代格式初值的选择对于迭代的收敛速度是至关重要的，当初值接近真实值的时候，迭代次数减少。
●在x0∈（-1，δ-）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：x。
计算结果显示，迭代序列局部收敛于1.730251，即根*3●在x0∈（δ-，δ）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：x。
由newton1.m的运行过程表明，在整个区间上均收敛于0，即根*2●在x0∈（δ，1）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：x。
计算结果显示，迭代序列局部收敛于-1.732051，即根*1在x0∈（1，+∞）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：x。
结果显示，以上初值迭代序列均收敛于1.732051，即根*3综上所述：(-∞，-1)区间收敛于-1.73205，(-1，δ)区间局部收敛于1.73205，局部收敛于-1.73205，(-δ，δ)区间收敛于0，(δ，1)区间类似于(-1，δ)区间，(1,∞)收敛于1.73205。
3、心得牛顿迭代法对于初值的选择要求较高，因此，在牛顿迭代时可现通过简单迭代法寻找相对准确一些的值来进行牛顿迭代。
对于方程有多解的问题，Newton法求方程根时，牛顿迭代要考虑局部收敛的问题，迭代序列收敛于某一个根有一定的区间限制，在一个区间上，可能会局部收敛于不同的根。
作业3、列主元素Gauss消去法对于某电路的分析，归结为求解线性方程组RI=V。
32 -13 0 0 0 -10 0 0 0-13 35 -9 0 -11 0 0 0 00 -9 31 -10 0 0 0 0 0R= 0 0 0 -30 57 -7 0 -5 00 0 0 0 -7 47 -30 0 00 0 0 0 0 -30 41 0 00 0 0 0 -5 0 0 27 -20 0 0 -9 0 0 0 -2 29V T=[-15，27，-23，0，-20,12，-7,7,10]T（1）编制解n阶线性方程组Ax=b的列主元Gauss消去法的通用程序；（2）用所编程序解线性方程组RI=V，并打印出解向量，保留5位有效数字；（3）在本编程之中，你提高了那些编程能力。
1、列主元Gauss消去法的通用程序函数：找每列的主元的函数function B=zhuyuan(B,t,N,M)for i=0:N-1-tif B(N-i,t)>B(N-i-1,t)c=zeros(1,M);for j=1:Mc(j)=B(N-i,j);B(N-i,j)=B(N-i-1,j);B(N-i-1,j)=c(j);endendend进行列消去的函数function B=xiaoqu(B,t,N,M)for i=t+1:Nl=B(i,t)/B(t,t);for j=t:MB(i,j)=B(i,j)-l*B(t,j);endend进行三角矩阵下的解函数function X=jie(X,B,N,M)for i=1:N-1s=B(N-i,M);for j=N-i+1:Ns=s-B(N-i,j)*X(j);endX(N-i)=s/B(N-i,N-i);end执行主程序：N=input('请输入线性方程组的阶数: N=');M=input('请输入增广矩阵阶数: M=');b=zeros(1,N);A=zeros(N,N);A=input('请输入系数矩阵:');b(1,:)=input('请输入方程组的右端向量:');b=b';B=[A,b];for t=1:N-1 B=zhuyuan(B,t,N,M);B=xiaoqu(B,t,N,M);endX=zeros(N,1);X(N)=B(N,M)/B(N,N);X=jie(X,B,N,M);X2、解题中线性方程组执行程序，输入矩阵A（即题中的矩阵R）和列向量b(即题中的V),如下：请输入线性方程组的阶数: n=9请输入增广矩阵阶数: m=10请输入系数矩阵A：A=[31,-13,0,0,0,-10,0,0,0,-15;-13,35,-9,0,-11,0,0,0,0,27;0,-9,31,-10,0,0,0,0,0,-23;0,0,-10,79,-30,0,0,0,-9,0;0,0,0,-30,57,-7,0,-5,0,-20;0,0,0,0,-7,47,-30,0,0,12;0,0,0,0,0,-30,41,0,0,-7;0,0,0,0,-5,0,0,27,-2,7;0,0,0,-9,0,0,0,-2,29,10];请输入方程组的右端向量b：[-15 27 -23 0 -20 12 -7 7 10]得到如下结果：A =31 -13 0 0 0 -10 0 0 0 -15-13 35 -9 0 -11 0 0 0 0 270 -9 31 -10 0 0 0 0 0 -230 0 -10 79 -30 0 0 0 -9 00 0 0 -30 57 -7 0 -5 0 -200 0 0 0 -7 47 -30 0 0 120 0 0 0 0 -30 41 0 0 -70 0 0 0 -5 0 0 27 -2 70 0 0 -9 0 0 0 -2 29 10x =-0.289230.34544-0.71281-0.22061-0.430400.15431-0.0578230.201050.290233、心得列主元Gauss最重要的就是如何通过找到最大主元，并交换行，如何进行消去，这需要很细心的循环，和很复杂的嵌套，通过本题的编程，感觉个人对于这种多层嵌套循环的处理能力提高了很多。