数值分析上机作业
数值分析上机作业(总)
数值分析上机实验一、解线性方程组直接法（教材49页14题）追赶法程序如下：function x=followup(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主程序如下：function zhunganfaA=[2 -2 0 0 0 0 0 0;-2 5 -2 0 0 0 0 0;0 -2 5 -2 0 0 0 0;0 0 -2 5 -2 0 0 0;0 0 0 -2 5 -2 0 0;0 0 0 0 -2 5 -2 0;0 0 0 0 0 -2 5 -2;0 0 0 0 0 0 -2 5];b=[220/27;0;0;0;0;0;0;0];x=followup(A,b)计算结果：x =8.14784.07372.03651.01750.50730.25060.11940.0477二、解线性方程组直接法（教材49页15题）程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:nfor i=1:1:nA(i,j)=(1+0.1*i)^(j-1);endendc=cond(A)d=rcond(A)当n=5时c =5.3615e+005d =9.4327e-007当n=10时c =8.6823e+011d =5.0894e-013当n=20时c =3.4205e+022d =8.1226e-024备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。
数值分析上机作业(MATLAB)
将系数矩阵 A 分解为：A=L+U+D
Ax=b
⇔ (D + L +U)x = b ⇔ Dx = −(L + U )x + b ⇔ x = −D −1(L + U )x + D −1b x(k +1) = −D −1 (L + U ) x(k ) + D −1b
输入 A，b 和初始向量 x
迭代矩阵 BJ , BG
否
ρ(B) < 1?
按雅各比方法进行迭代
否
|| x (k+1) − x(k) ||< ε ?
按高斯-塞德尔法进行迭代
否
|| x(k+1) − x (k ) ||< ε ?
输出迭代结果
图 1 雅各布和高斯-赛德尔算法程序流程图
1.2 问题求解
按图 1 所示的程序流程，用 MATLAB 编写程序代码，具体见附录 1。解上述三个问题 如下
16
-0.72723528355328
0.80813484897616
0.25249261987171
17
-0.72729617968010
0.80805513082418
0.25253982509100
18
-0.72726173942623
0.80809395746552
0.25251408253388
0.80756312717373
8
-0.72715363032573
0.80789064377799
9
-0.72718652854079
数值分析大作业
数值分析上机作业（一）一、算法的设计方案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 。
数值分析上机作业
第一章第二题(1) 截断误差为104-时：k=1;n=0;m=0;x=0;e=1e-4;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141792613595791m =5001（2）截断误差为108-时：k=1;n=0;m=0;x=0;e=1e-8;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141592673590250m =50000001由以上计算可知，截断误差小于104-时，应取5001项求和，π=3.141792613595791；截断误差小于108-时，应取50000001项求和，π=3.141592673590250。
第二章第二题a=[0 -2 -2 -2 -2 -2 -2 -2];b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0];v=220;r=27;d=[v/r 0 0 0 0 0 0 0];n=8;for i=2:na(i)=a(i)/b(i-1);b(i)=b(i)-c(n-1)*a(i);d(i)=d(i)-a(i)*d(i-1);end;d(n)=d(n)/b(n);for i=n-1:-1:1d(i)=(d(i)-c(i)*d(i+1));end;I=d'I =1.0e+002 *1.490717294184090.704617906351300.311568212434910.128623612390290.049496991380330.017168822994210.004772412363470.00047741598598第三章第一题(1)Jacobi迭代法：b=[12;-27;14;-17;12]x = [0;0;0;0;0;]k = 0;r = 1;e=0.000001A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] D = diag(diag(A));B = inv(D)*(D-A);f = inv(D)*b;p = max(abs(eig(B)));if p >= 1'迭代法不收敛'returnendwhile r >ex0 = x;x = B*x0 + f;k = k + 1;r = norm (x-x0,inf);endxk计算结果：x =1.0000-2.00003.0000-2.00001.0000k =65(2) 高斯赛德尔迭代：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;]x=[0;0;0;0;0];b=[12;-27;14;-17;12]c=0.000001L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-L)*U*x+inv(D-L)*b;k=1;while norm(X-x,inf)>= cx=X;X=inv(D-L)*U*x+inv(D-L)*b;k=k+1;endXk计算结果：X =1.0000-2.00003.0000-2.00001.0000k =37(3) SOR：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15] x=[0;0;0;0;0];b=[12;-27;14;-17;12]e=0.000001w=1.44;L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*bn=1;while norm(X-x,inf)>=ex=X;X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b;n=n+1;endXn计算结果：X =1.0000-2.00003.0000-2.00001.0000n =22由以上可知，共轭梯度法收敛速度明显快于Jacobi法和G-S法。
东南大学出版社第二版《数值分析》上机作业答案(前三章)
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;
数值分析上机作业
第一章（1）按从大到小的顺序计算S的程序为：N#include<stdio.h>void main(){float n,j,s=0;printf("请输入n：");scanf("%f",&n);for(j=2;j<=n;j++)s=s+1/(j*j-1);printf("s=%f\n",s);}(2)按从小到大的顺序计算S的程序为：N#include<stdio.h>void main(){float n,j,s=0;printf("请输入n：");scanf("%f",&n);for(j=n;j>=2;j--)s=s+1/(j*j-1);printf("s=%f\n",s);}(3)精确值，从大到小，从小到大的计算程序为：#include<stdio.h>void main(){float n,j,s1,s2=0,s3=0;printf("请输入n：");scanf("%f",&n);s1=(3.0/2-1/n-1/(n+1))/2;for(j=2;j<=n;j++)s2=s2+1/(j*j-1);for(j=n;j>=2;j--)s3=s3+1/(j*j-1);printf("精确算法的结果为：%f\n",s1);printf("从大到小顺序的结果为：%f\n",s2); printf("从小到大顺序的结果为：%f\n",s3); }结果如下：精确值从大到小的值从小到大的值有效位数从大到小从小到大210S0.740049 0.740049 0.74005 6 5410S0.7499 0.749852 0.7499 4 4610S0.749999 0.749852 0.749999 3 6(4)通过本上机题，我们可以看出：当n较小时，两种算法的结果都很接近精确值，而当n较大时，两种算法的结果差别较大，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。
高等数值分析上机作业
画图分析
4
第 10 章 常微分方程的初边值问题解法
上机作业 2：
1. 分别用 Euler 方法（h=0.025）、改进的 Euler 方法（h=0.05）、经
典的 R-K 方法（h=0.1），求初值问题的数值计算 0.1,0.2,0.3,0.4,0.5
处解的近似值。
y y x2 1 y(0) 0.5
主程序： >> fun=inline('y-x.^2+1'); >> [x,y]=Euler1(fun,0.5,0.1)
运行结果：
x=
0.1000 0.2000 0.3000 0.4000 0.5000
y=
0.5000 0.6540 0.8200 0.9970 1.1841
用经典的R-K方法（）求解一阶微分方程组的初值问题，并分别画出
3
h=0.01;x=-1:h:1; B(1,1)=quad('asin(x)./(sqrt(1-x.^2))',-1,1); B(2,1)=quad('asin(x).*x./sqrt(1-x.^2)',-1,1); B(3,1)=quad('asin(x).*(2*x.^2-1)./(sqrt(1-x.^2))',-1,1); B(4,1)=quad('asin(x).*(4*x.^3-3*x)./(sqrt(1-x.^2))',-1,1); B A=diag([pi,pi/2,pi/2,pi/2]); c=inv(A)*B y=asin(x); plot(x,y,'r') hold on z=0.9204*x+0.4296*x.^3; plot(x,z,'k') hold on s=1.2733*x+0.1415*(4*x.^3-3*x); plot(x,s,'c') legend('y=arcsin(x)','y=0.9204x+0.4296x^3','s=1.2733x+0.141 5(4x^3-3*)') 最后可以得到切比雪夫拟合多项式为 y 1.2733x 0.1415(4x3 3x) ，
数值分析上机作业
《数值分析》上机作业（第一二三章）学院：电气工程学院班级：电气13级硕士2班教师：石佩虎老师姓名：**学号： ******第一章实验1 舍入误差与有效数设2211NN j S j==-∑，其精确值为1311()221N N --+。
(1) 编制按从大到小的顺序222111 (21311)N S N =+++---，计算N S 的通用程序； (2) 编制按从小到大的顺序222111...1(1)121N S N N =+++----，计算N S 的通用程序； (3) 按两种顺序分别计算210S 、410S 、610S ，并指出有效位数（编制程序时用单精度）； (4) 通过本上机题你明白了什么？解答如下：(1). 按从大到小的顺序计算N S 的通用程序如下所示： n=input('Please Input an N (N>1):'); y=0;accurate=1/2*(3/2-1/n-1/(n+1)); %精确值 for i=2:1:n %从大到小的顺序 x=1/(i^2-1);x=single(x); y=y+x; enderror= accurate-y; format long;disp('____________________________________________________'); disp('The value of Sn from large to small is:'); disp(y);disp('The value of error is:'); disp(error);(2) 编制按从小到大的顺序计算N S 的通用程序如下所示： n=input('Please Input an N (N>1):'); y=0;accurate=1/2*(3/2-1/n-1/(n+1)); for i=n:-1:2 x=1/(i^2-1);x=single(x); y=y+x;enderror= accurate-y; format long;disp('____________________________________________________'); disp('The value of Sn from large to small is:'); disp(y);disp('The value of error is:'); disp(error);(3) 计算结果：按从大到小的顺序计算得：(4)总结：当我们采用不同的计算顺序，对于同一个计算式，会得出不同的结果。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第二次上机作业
一. 任务：
用MATLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。
并用此
程序进行数值试验，写出计算实习报告。
二. 程序功能要求：
在后面的附一leastp.m 的基础上进行修改，使其更加完善。
要求算法程序可以适应不同的具体函数，具有一定的通用性。
所编程序具有以下功能：
1. 用Lengendre 多项式做基，并适合于构造任意次数的最佳平方逼近多项式。
可利用递推关系 0112()1,()()(21)()(1)()/2,3,.....
n n n P x P x x
P x n xP x n P x n n --===---⎡⎤⎣⎦=
2. 被逼近函数f(x)不用内联函数构造，而改用M 文件建立数学函数。
这样，此程序可通过修
改建立数学函数的M 文件以适用不同的被逼近函数（要学会用函数句柄）。
3. 要考虑一般的情况]1,1[],[)(+-≠∈b a x f 。
因此，程序中要有变量代换的功能。
4. 计算组合系数时，计算函数的积分采用变步长复化梯形求积法（见附三）。
5. 程序中应包括帮助文本和必要的注释语句。
另外，程序中也要有必要的反馈信息。
6. 程序输入：（1）待求的被逼近函数值的数据点0x (可以是一个数值或向量)
（2）区间端点：a,b 。
7. 程序输出：（1）拟合系数：012,,,...,n c c c c
（2）待求的被逼近函数值
00001102200()()()()()n n s x c P x c P x c P x c P x =++++
三：数值试验要求：
1. 试验函数：()cos ,[0,4]f x x x x =∈+；也可自选其它的试验函数。
2. 用所编程序直接进行计算，检测程序的正确性，并理解算法。
3. 分别求二次、三次、。
最佳平方逼近函数）x s (。
4. 分别作出逼近函数）x s (和被逼近函数)(x f 的曲线图进行比较。
（分别用绘图函数plot(0x ,s(0x ))和fplot(‘x cos x ’,[x 1 x 2,y 1,y 2])） 四：计算实习报告要求：
1.简述方法的基本原理，程序功能，使用说明。
2.程序中要加注释。
3.对程序中的主要变量给出说明。
4.附源程序及计算结果。
附： 一、 参考程序
Lengendre 多项式作基的函数最佳平方逼近算法程序LEASTP.m
（此程序只适用于对函数x
xe x f =)( 构造最佳平方逼近多项式）
function [c,s]=leastp(x)
%LEASTP.m:least-square fitting with legendre polynomials
p1=1;
p2=inline('x','x'); p3=inline('(3*x^2-1)','x'); pp1=1;
pp2=inline('x.^2','x');
pp3=inline( '(3*x.^2-1)/2.*(3*x.^2-1)/2','x'); fp1=inline('x.*exp(x)','x'); fp2=inline('(x.^2).*exp(x)','x'); fp3=inline('(x.*exp(x)).*(3*x.^2-1)/2','x'); c(1)=quad(fp1,-1,1)/2;
c(2)=quad(fp2,-1,1)/quad(pp2,-1,1); c(3)=quad(fp3,-1,1)/quad(pp3,-1,1); s=c(1)+c(2)*p2(x)+c(3)*p3(x);
二、被逼近函数用M 文件建立（例如下面）
function f=fun(x) f=1./(1+x.^2);
三、变步长复化梯形求积公式的算法
()11111.,(()())22.0,2
3.(),3.15.*2
6.,
7.,,2
b a
h b a T f a f b h
H x a H H f x x x h
x b T T h H T T I T I h
h T T ε-=-=
+==+
=+=+<=+-<==
=4.若，则转若则，输出,停机。
转2. 121
0[()()]2
11
222((21))221,2,n n n n n j b a
T f a f b H T T T b a b a
f a j n n n -=-=
+=+=+
--++=∑
第三次上机作业（题目自拟）结合自己专业选择适当的数值方法和数值算例（注重实际应用）上机作业要求
•简述方法的基本原理。
•程序功能，变量说明，使用说明。
•程序中要加语句注释。
•自选数值例子，检验程序的正确性。
•附原程序及计算结果复印件(A4纸)。
•将上机报告按顺序装订成册。