数值分析第二章上机
数值分析上机作业(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
数值分析---上机考试大纲
东华大学研究生《数值分析》实验考试大纲教材：«数值分析及其MA TLAB实验»姜健飞胡良剑唐俭编考试规则领座试卷不同，开卷，解答全部用笔写在考卷上，作图题只需手画草图。
开考前可将准备程序Copy到硬盘, 但是开考后不允许用软盘，也不允许上网。
评分原则满分20。
对一题得6分, 对两题得11分, 对三题得15分. 对四题得18分. 基本正确题酌情给分。
类型1：使用Matlab命令的计算题共3题主要使用如下MA TLAB命令解题：第二章（1）用矩阵除法解线性方程组；（2）行列式det、逆inv；（3）特征值、特征向量eig；（4）范数和条件数；第三章（1）用roots求多项式的根；（2）用fzero解非线性方程；（3）用fsolve解非线性方程组；第四章（1）多项式插值和拟合polyfit(2) 线性插值interp1(3) 样条插值spline, csape（4）最小二乘拟合lsqcurvefit第五章（1）用diff或gradiet求导数（2）用trapz、quad或quadl求积分；（3）用dblquad或triplequad求重积分；第六章（1）用ode45求解微分方程；（2）用ode45求解微分方程组；（3）用ode45求解高阶微分方程；类型2：使用课本程序的计算题共1题（不必将课本程序部分写在考卷上）第二章nagauss nagauss2 nalu nalupad第三章nabisect nanewton nags naspgs nasor第四章nalagr naspline nafit naorthfit第五章natrapz nagsint naromberg naadapt dblquad2第六章naeuler naeulerb naeuler2 nark4 nark4v naeuler2s类型3：编程题共1题(必须将程序写在考卷上）要求使用MA TLAB控制流语句编程，主要涉及for, while, if等语句以及关系与逻辑运算，M 函数编写。
东南大学数值分析上机题答案
东南⼤学数值分析上机题答案数值分析上机题第⼀章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）感想：通过本上机题，我明⽩了，从⼩到⼤计算数值的精确位数⽐较⾼⽽且与真值较为接近，⽽从⼤到⼩计算数值的精确位数⽐较低。
东南大学数值分析上机报告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 绪论 (1)2 实验题目(一) (2)2.1 题目要求 (2)2.2 NEWTON插值多项式 (3)2.3 数据分析 (4)2.3.1 NEWTON插值多项式数据分析 (4)2.3.2 NEWTON插值多项式数据分析 (6)2.4 问答题 (6)2.5 总结 (7)3 实验题目(二) (8)3.1 题目要求 (8)3.2 高斯-塞德尔迭代法 (8)3.3 高斯-塞德尔改进法—松弛法 (9)3.4 松弛法的程序设计与分析 (9)3.4.1 算法实现 (9)3.4.2 运算结果 (9)3.4.3 数据分析 (11)4 实验题目(三) (13)4.1 题目要求 (13)4.2 RUNGE-KUTTA 4阶算法 (13)4.3 RUNGE-KUTTA 4阶算法运算结果及数值分析 (14)总结 (16)附录A (17)1绪论数值分析是计算数学的一个主要部分，它主要研究各类数学问题的数值解法，以及分析所用数值解法在理论上的合理性。
实际工程中的数学问题非常复杂，所以往往需要借助计算机进行计算。
运用数值分析解决问题的过程：分析实际问题，构建数学模型，运用数值计算方法，进行程序设计，最后上机计算求出结果。
数值分析这门学科具有面向计算机、可靠的理论分析、好的计算复杂性、数值实验、对算法进行误差分析等特点。
本学期开设了数值分析课程，该课程讲授了数值分析绪论、非线性方程的求解、线性方程组的直接接法、线性方程组的迭代法、插值法、函数逼近与曲线拟合、数值积分和数值微分、常微分方程初值问题的数值解法等内容。
其为我们解决实际数学问题提供了理论基础，同时我们也发现课程中很多问题的求解必须借助计算机运算，人工计算量太大甚至无法操作。
所以学好数值分析的关键是要加强上机操作，即利用计算机程序语言实现数值分析的算法。
本报告就是基于此目的完成的。
本上机实验是通过用计算机来解答数值分析问题的过程，所用的计算工具是比较成熟的数学软件MATLAB。
《数值分析》第二章上机习题
上机实习要求: 编程可以用C、C++、Matlab, 但不允许使用内置函数完成主要功能。
第2 章1. 已知函数表如下x 10 11 12 13ln(x) 2.3026 2.3979 2.4849 2.5649试分别用线性插值与二次插值计算ln11.75 的近似值，并估计截断误差.2. 已知函数表x 0.1 0.2 0.3 0.4sin(x) 0.09983 0.19867 0.29552 0.38942试分别用Newton 前插与后插公式（1、2、3 阶）计算sin(0.22)的近似值。
要求，比较所得结果，思考如何选取节点。
3．构造函数表cos(x)：已知节点x k=k⋅/20 (k=0,1,…,20)处的函数值. 用一次和二次Lagrange 插值公式求cos(x)在x k_i(i=1,2,3)( x k_i=x k+( x k +1- x k)/4⋅i). 请用你计算的值连成函数图形，与标准图形比较。
4．已知直升飞机旋转机翼外形曲线轮廓上的某些型值点（见表），及端点处的一阶导数值y'(x) =. , y'(x) = .0 1 86548 18 0 046115试计算该曲线上横坐标为2,4,6,12,16,30,60,110,180,280,400,515 处的纵坐标（要求该曲线具有二阶光滑度）.k 0 1 2 3 4 5 6x k 0.52 3.1 8.0 17.95 28.65 39.62 50.65y k 5.28794 9.4 13.84 20.2 24.9 28.44 31.1k 7 8 9 10 11 12 13x k 78 104.6 156.6 208.6 260.7 312.5 364.4y k 35 36.5 36.6 34.6 31.0 26.34 20.9k 14 15 16 17 18x k 416.3 468 494 507 520y k 14.8 7.8 3.7 1.5 0.2。
数值分析上机实验指导书
“数值计算方法”上机实验指导书实验一 误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。
对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。
通过本实验可获得一个初步体会。
数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。
病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。
问题提出：考虑一个高次的代数多项式)1.1()()20()2)(1()(201∏=−=−−−=k k x x x x x p显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。
现考虑该多项式的一个扰动)2.1(0)(19=+x x p ε其中ε是一个非常小的数。
这相当于是对（1.1）中19x 的系数作一个小的扰动。
我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。
实验内容：为了实现方便，我们先介绍两个MATLAB 函数：“roots ”和“poly ”。
roots(a)u =其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。
设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程01121=+++++−n n n n a x a x a x a的全部根；而函数 poly(v)b =的输出b 是一个n+1维向量，它是以n 维向量v 的各分量为根的多项式的系数（从高到低排列）。
可见“roots ”和“poly ”是两个互逆的运算函数。
))20:1((;)2();21,1(;000000001.0ve poly roots ess ve zeros ve ess +===上述简单的MATLAB 程序便得到（1.2）的全部根，程序中的“ess ”即是（1.2）中的ε。
数值分析上机作业
《数值分析》上机作业（第一二三章）学院：电气工程学院班级：电气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)总结：当我们采用不同的计算顺序，对于同一个计算式，会得出不同的结果。
数值分析上机实验报告
数值分析上机实验报告导言：本次上机实验主要是针对数值分析课程中的一些基本算法进行实验验证。
实验内容包括迭代法、插值法、数值积分和常微分方程的数值解等。
在实验过程中，我们将会使用MATLAB进行算法的实现，并对结果进行分析。
一、迭代法迭代法是解决函数零点、方程解等问题的常用方法。
我们将选择几个常见的函数进行迭代求根的实验。
（1）二分法二分法是一种简单而有效的迭代求根法。
通过函数在区间两个端点处的函数值异号来确定函数在区间内存在零点，并通过不断缩小区间来逼近零点。
（2）牛顿法牛顿法利用函数的一阶导数和二阶导数的信息来逼近零点。
通过不断迭代更新逼近值，可以较快地求得零点。
实验结果表明，对于简单的函数，这两种迭代法都具有很好的收敛性和稳定性。
但对于一些复杂的函数，可能会出现迭代失效或者收敛速度很慢的情况。
二、插值法插值法是在给定一些离散数据点的情况下，通过构造一个插值函数来逼近未知函数的值。
本实验我们将使用拉格朗日插值和牛顿插值两种方法进行实验。
（1）拉格朗日插值拉格朗日插值通过构造一个多项式函数来逼近未知函数的值。
该多项式经过离散数据点，并且是唯一的。
该方法简单易懂，但插值点越多，多项式次数越高，插值函数的精度也就越高。
（2）牛顿插值牛顿插值利用差商的概念，通过构造一个插值多项式来逼近未知函数的值。
与拉格朗日插值相比，牛顿插值的计算过程更加高效。
但同样要求插值点的选择要合理，否则可能出现插值函数不收敛的情况。
实验结果表明，这两种插值方法都能够很好地逼近未知函数的值。
插值点的选择对插值结果有很大的影响，过多或者过少的插值点都可能导致插值结果偏离真实函数的值。
三、数值积分数值积分是一种将定积分问题转化为数值求和的方法。
本实验我们将使用复合梯形求积法和复合辛普森求积法进行实验。
（1）复合梯形求积法复合梯形求积法将定积分区间等分为若干小区间，然后使用梯形公式对每个小区间进行近似求积，最后将结果相加得到整个定积分的近似值。
《应用数值分析》第二次上机作业
可利用递推关系
P ( x ) 1, P ( x ) x
0
1
P( n
x
) ( 2 n1)xP ( n1
x)
(n
1)P ( n2
x ) /
n
n 2, 3, .....
2、被逼近函数 f(x)不用内联函数构造，而改用 M 文件建立数学函数。这样，
此程序可通过修改建立数学函数的 M 文件以适用不同的被逼近函数（要
fa=feval(func,a);fb=feval(func,b); T2n=(fa+fb)*(b-a)/2.0; else n=2^(k-2); h=(b-a)/n; x=a+h/2; sum=0.0; for i=1:n
fx=feval(func,x); sum=sum+fx;
x=x+h;
end
二、被逼近函数用 M 文件建立（例如下面） function f=fun(x) f=1./(1+x.^2);
三、变步长复化梯形求积公式的算法
ba
1.h b a ,T
( f (a ) f (b ))
2
h 2.H 0, x a
2
3.H H f ( x ), x x h
2、三次最佳平方逼近函数：
输入 syms x;
[c,s]=leastp_ww(x,6,3,-1,1)
运行结果：
s=
17336850270911903/18014398509481984-9565132890499617/18014398509481984
*x^2
c=
0.785398004451730 0 -0.353981025621068
《数值分析》上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
数值分析上机作业(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,用试算法确定最佳松弛因子。
《数值分析》上机实验报告
数值分析上机实验报告x k x k - f(X k) f (X k)《数值分析》上机实验报告1. 用Newt on法求方程X7-X4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001 ）。
1.1理论依据：设函数在有限区间［a，b］上二阶导数存在，且满足条件1. f(x)f(b) 02. f(x)在区间［a, b］上不变号3f(x) = 0;4」f (c)〔f .(x) |,其中c是a,b中使mir(| f .(a), f .(b) |)达到的一个b -a则对任意初始近似值x0• ［a,b］,由Newton迭代过程込f(x k )X“ M(Xk) = Xk — T^,k = 0,1,2,3…f'(X k)所生的迭代序列 % ［平方收敛于方程f(x)=0在区间［a,b］上的惟一解: 令7 4f(x)=x -28x 14, f (0.1) 0, f(1.9) ：：0f (x) =7x6-112x3=7x3(x3-16) ：：： 0f (x) =42x5-336x2=42x2(x3-8) ：： 0f (1.9) f (1.9) 0故以1.9为起点x0 =1.9如此一次一次的迭代，逼近X的真实根。
当前后两个的差<=出寸，就认为求出了近似的根。
本程序用Newton法求代数方程（最高次数不大于10）在（a,b ）区间的根//限制循环次数1.2 C 语言程序原代码：#i nclude<stdio.h> #in clude<math.h> mai n() {double x2,f,f1; double x1=1.9; // 取初值为 1.9do {x2=x1;f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3); x 仁 x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); printf("计算结果：x=%f\n",x1);}1.3运行结果:* D:\VC + +\EXERCIS E\Debu g\l1.4 MATLAB上机程序fun cti on y=Newt on( f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2; breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))v=epsd=1; breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey=奇异’endfun cti on y=df(x)y=7*x A6-28*4*x A3;Endfunction y=f(x) y=x A7-28*x A4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newto n('f,'df,x0,eps,M);>> vpa(x,7)1.5问题讨论：1•使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
数值分析上机实验报告
数值分析上机实验报告数值分析上机实验报告一、引言数值分析是一门研究利用计算机进行数值计算的学科。
通过数值分析，我们可以使用数学方法和算法来解决实际问题，例如求解方程、插值和逼近、数值积分等。
本次上机实验旨在通过编程实现数值计算方法，并应用于实际问题中。
二、实验目的本次实验的目的是掌握数值计算方法的基本原理和实现过程，加深对数值分析理论的理解，并通过实际应用提高编程能力。
三、实验内容1. 数值求解方程首先，我们使用二分法和牛顿迭代法分别求解非线性方程的根。
通过编写程序，输入方程的初始值和精度要求，计算得到方程的根，并与理论解进行对比。
2. 数值插值和逼近接下来，我们使用拉格朗日插值和最小二乘法进行数据的插值和逼近。
通过编写程序，输入给定的数据点，计算得到插值多项式和逼近多项式，并绘制出插值曲线和逼近曲线。
3. 数值积分然后，我们使用梯形法和辛普森法进行定积分的数值计算。
通过编写程序，输入被积函数和积分区间，计算得到定积分的近似值，并与解析解进行比较。
四、实验步骤1. 数值求解方程（1）使用二分法求解非线性方程的根。
根据二分法的原理，编写程序实现二分法求解方程的根。
（2）使用牛顿迭代法求解非线性方程的根。
根据牛顿迭代法的原理，编写程序实现牛顿迭代法求解方程的根。
2. 数值插值和逼近（1）使用拉格朗日插值法进行数据的插值。
根据拉格朗日插值法的原理，编写程序实现数据的插值。
（2）使用最小二乘法进行数据的逼近。
根据最小二乘法的原理，编写程序实现数据的逼近。
3. 数值积分（1）使用梯形法进行定积分的数值计算。
根据梯形法的原理，编写程序实现定积分的数值计算。
（2）使用辛普森法进行定积分的数值计算。
根据辛普森法的原理，编写程序实现定积分的数值计算。
五、实验结果与分析1. 数值求解方程通过二分法和牛顿迭代法，我们成功求解了给定非线性方程的根，并与理论解进行了对比。
结果表明，二分法和牛顿迭代法都能够较好地求解非线性方程的根，但在不同的问题中，二者的收敛速度和精度可能会有所差异。
数值分析上机实习题
数值分析上机实习题第2章插值法1. 已知函数在下列各点的值为试⽤四次⽜顿插值多项式）（x p 4及三次样条韩式)(S x (⾃然边界条件)对数据进⾏插值。
⽤图给出(){}10,11,1,0,08.02.0,,x i =+=i x y i i ，）（x p 4及)(x S Python 代码import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.font_manager import FontPropertiesfont_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=12) #求⽜顿n 次均差 def qiujuncha(x,f,n): for i in range(1,n): for j in range(4,i-1,-1):f[j]= (f[j] - f[j-1])/(x[j]-x[j-i]) #根据⽜顿多项式求值 def niudun(x,f,x1): sum = f[0]; tmp = 1;for i in range(1,5): tmp *= (x1-x[i-1]) sum = sum + f[i]*tmp return sum#⽜顿插值画图 def drawPic(x,f):x1 = np.linspace(0.2, 1, 100) plt.plot(x1, niudun(x,f,x1))plt.title(u"⽜顿四次插值",fontproperties=font_set) plt.xlabel(u"x 轴",fontproperties=font_set) plt.ylabel(u"y 轴", fontproperties=font_set) plt.show() def qiu_h(x,h): n = len(x) -1 for i in range(n): print(i)h[i] = x[i+1]-x[i]#⾃然边界条件下的三次样条插值求Mdef qiu_m(h,f,o,u,d):n = len(h)o[0] = 0u[n] = 0d[n] = d[0] = 0a = []for i in range(1,n):u[i] = h[i-1]/(h[i-1]+h[i])for i in range(1,n):o[i] = h[i]/(h[i-1]+h[i])for i in range(1,n-1):d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])t = [0 for i in range(5)]t[0] =2t[1] = o[0]a.append(t)for i in range(1,n):t = [0 for i in range(5)]t[i - 1] = u [i + 1]t[i] = 2t[i + 1] = o [i + 1]a.append(t)t = [0 for i in range(5)]t[n - 1] = u[n]t[n] = 2a.append(t)tmp = np.linalg.solve(np.array(a),np.array(d))m = []for i in range(5):m.append(tmp[i])return m#根据三次条插值函数求值def yangtiao(x1,m,x,y,h,j):returnm[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] def main():x = [0.2, 0.4, 0.6, 0.8, 1.0]y = [0.98, 0.92, 0.81, 0.64, 0.38]f = y[:]f1 = y[:]h = [0.2,0.2,0.2,0.2]u = [0 for n in range(5)]d = [0 for n in range(5)]o = [0 for n in range(5)] qiujuncha(x,f,4) qiujuncha(x,f1,2)m = qiu_m(h,f1,o,u,d) x1 = np.linspace(0.2, 0.4, 10)p1= plt.plot(x1, yangtiao(x1,m,x,y,h,0),color='red') x1 = np.linspace(0.4, 0.6, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 1),color='red') x1 = np.linspace(0.6, 0.8, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 2),color='red') x1 = np.linspace(0.8, 1.0, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 3),color='red') x1 = np.linspace(0.2, 1.0, 40)p2 = plt.plot(x1,niudun(x,f,x1),color='green') plt.xlabel(u"x 轴", fontproperties=font_set) plt.ylabel(u"y 轴",fontproperties=font_set) plt.title("三次样条插值和⽜顿插值")plt.legend(labels=[u'三次样条插值',u'⽜顿插值'],prop=font_set,loc="best") plt.show() main()实验结果运⾏结果可得插值函数图（如图1-1），4次⽜顿插值函数）（x p 4和三次样条插值函数)(x S 如下：)6.0(*)4.0(*)2.0(625.0)4.0(*)2.0(*3.098.0)(4-------=x x x x x x x P 98.0)8.0(*)6.0(*)4.0(*)2.0(*20833.0+-----x x x x]4.0,2.0[),2.0(467.4)4.0(9.4)2.0(167.1)(S 3∈-+-+-=x x x x x]6.0,4.0[),4.0(113.4)6.0(6467.4)4.0(575.1)6.0(167.1)(S 33∈-+-+----=x x x x x x ]8.0,6.0[),6.0(2.3)8.0(113.4)6.0(575.1)(S 3∈-+-+--=x x x x x]0.1,8.0[),8.0(9.1)0.1(2.3)(S ∈-+-=x x x x图1-1三次样条插值和⽜顿插值图2.在区间[-1,1]上分别取n = 10,20⽤两组等距节点对龙格函数做多项式插值三次样条插值，对每个n值画出插值函数及图形。
东南大学-数值分析-第二章-牛顿迭代法
东南大学-数值分析-第二章-牛顿迭代法第二章非线性方程的解法某某某某（学号）某某某某（姓名）算法与程序题目见教材P56上机题目20。
一、算法原理根据题目的要求，是关于用牛顿迭代法法求解方程f(某)0的通用算法。
该法是一种通过斜率迭代的算法，其速度比二分法和简单迭代法都要快。
其简单原理如下：设fC2[a,b],且存在数p[a,b],满足f(p)0。
如果f(p)0,则存在一个数0,对任意初始值p0[p,p],使得由如下定义的迭代序列{pk}k0收敛到p：pkg(pk1)pk1f(pk1),其中k1,2,f(pk1)(1)对于函数f(某)某3/3某=0，则其递推规则是32pkpk21,其中k1,2,3pk1-3(2)定义序列{pk}则序列{pk}也可表示为limpk某现简要证明：k0，k0收敛到某，某对于f(某)某3/3某，得f'(某)某2-1，写出牛顿迭代公式f(某)某3/3某g(某)某某2f(某)某-1(3)该公式可化简为2某3g(某)23某3(4)二、流程图题目要求于用牛顿迭代法法求解方程f(某)0的通用算法。
其计算过程主要第二章非线性方程的解法用到迭代g(某)某f(某)，图流程图1所示。
f(某)输入各参数k=1迭代pkg(pk1)pk1f(pk1),其中k1,2,f(pk1)Tbreak计算各误差误差在允许范围之内Fk=k+1k三、计算代码核心代码1)p1=……;2)if(err程序1：Newton.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Decription:牛顿迭代法%Author:panyunqiang%Veroin:1.0%Date:2022-9-21%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%f unction[p0,err,k,y]=Newton(p0,delta,epilon,ma某N)%input-p0itheinitialappro某imationtoazerooff%-deltaithetoleranceforp0%-epilonithetoleranceforthefunctionvaluey%-ma某Nithema某iumnumberofiteration%output-p0itheNewtonappro某imationtoazero%-erritheerroretimateforp0东南大学《数值分析》上机练习——算法与程序设计实验报告%-kithenumberofiteration%-yithefunctionvaluef(p0)fork=1:ma 某N%%递归p1=2某p0^3/(3某p0^2-3);%%计算误差err=ab(p1-p0);relerr=2某err/(ab(p1)+delta);p0=p1;%%当前求出的根的函数值y=p0^3/3-p0;%%判断if(err程序2：Newton_Step.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%Decription:寻找题目中关于牛顿迭代法收敛的尽可能大的delta%搜索步进为tep=10^(-6),即精确到小数点后六位%Author:panyunqiang%Veroin:1.0%Date:2022-9-21%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %formatlongtep=10^(-6);delta=10^-8;epilon=10^-8;ma某N=1000;p=0.6;[p0,err,k,y]=Newton(p,delta,epilon,ma某N);while((ab(p0)<=epilon)&(p0~=NaN))p=p+tep;[p0,err,k,y]=Newton(p,delta,epilon,ma某N);endp-tep四、计算结果及分析a)运行程序Newton_Step.m，获得Newton局部收敛于某2=0的初始值的范围=0.774596，六位有效数字。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
用列主元消元法解线性方程组AX=b的MATLAB程序function[RA,RB,n,X]=liezhu(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA~=RB，所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RB<n，所以此方程组有无穷多解.') endend习题3.33.在MATLAB工作窗口输入程序>> A=[1 1 1;1 3 9;1 7 49];b=[6;5;2];[RA,RB,n,X]=liezhu(A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.RA = 3RB = 3n = 3X = 6.3750-0.3333-0.0417将矩阵A进行直接LU分解的MATLAB程序function hl=zhjLU(A)[n n]=size(A);RA=rank(A);if RA~=ndisp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'),RA,hl=det(A); returnendif RA==nfor p=1:nh(p)=det(A(1:p,1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl 依次如下：'),hl;RAreturnendendif h(1,i)~=0disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend习题3.41.(1) 在MATLAB工作窗口输入程序>> A=[2 4 -6;1 5 3;1 3 2];hl=zhjLU(A)运行后输出结果请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 3U =2.0000 4.0000 -6.00000 5.0000 6.00000 0 3.8000L =1.0000 0 00.5000 1.0000 00.5000 0.2000 1.0000hl =2 6 18(2) 在MATLAB工作窗口输入程序>> A=[1 1 6;-1 2 9;1 -2 3];hl=zhjLU(A)运行后输出结果请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA =3U = 1.0000 1.0000 6.00000 2.0000 15.00000 0 19.5000L = 1.0000 0 0-1.0000 1.0000 01.0000 -1.5000 1.0000hl = 1 3 36用P范数讨论AX=b解和A的性态的MATLAB程序function Acp=zpjwc(A,jA,b,jb,p)Acp=cond(A,p);dA=det(A);X=A\b;dertaA=A-jA;PndA=norm(dertaA,p);dertab=b-jb;Pndb=norm(dertab,p);if Pndb>0jX=A\jb;Pnb=norm(b,p);PnjX=norm(jX,p);dertaX=X-jX;PnjdX=norm(dertaX,p);jxX=PnjdX/PnjX;PnjX=norm(jX,p);PnX=norm(X,p);jxX=PnjdX/PnjX;xX=PnjdX/PnX;Pndb=norm(dertab,p);xAb=Pndb/Pnb;Pnbj=norm(jb,p);xAbj=Pndb/Pnbj;Xgxx=Acp*xAb;endif PndA>0jX=jA\b;dertaX=X-jX;PnX=norm(X,p);PnjdX=norm(dertaX,p);PnjX=norm(jX,p);jxX=PnjdX/PnjX;xX=PnjdX/PnX;PnjA=norm(jA,p);PnA=norm(A,p);PndA=norm(dertaA,p);xAbj=PndA/PnjA;xAb=PndA/PnA;Xgxx=Acp*xAb;endif (Acp>50)&&(dA<0.1)disp('请注意：AX=b是病态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X,jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意： AX=b是良态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'end习题3.61.在MATLAB工作窗口输入程序>>A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];jA=A;b=[32 23 33 31]';jb=[32.1 22.9 22.2 30.9]';Acp=zpjwc(A,jA,b,jb,1) 运行后输出结果请注意：AX=b是良态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：Acp =4.4880e+03dA =1.0000ans =1.0000 1.0000 1.0000 1.0000ans =-99.8000 172.7000 -50.0000 31.6000xX =88.5250jxX =1.0000Xgxx =418.6286xAb =0.0933xAbj =0.1027Acp =4.4880e+03用雅可比迭代解线性方程组AX=b的MATLAB主程序function X=jacdd(A,b,X0,P,wucha,max1)[n m]=size(A);for j=1:ma(j)=sum(abs(A(:,j)))-2*(abs(A(j,j)));endfor i=1:nif a(i)>=0disp('请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛')returnendendif a(i)<0disp('请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛') endfor k=1:max1k;for j=1:mX(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1: j-1,j+1:m]))/A(j,j);endX;djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\b;if (djwcX<wucha)&&(xdwcX<wucha)disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：') breakendendif (djwcX>wucha)&&(xdwcX>wucha)disp('请注意：雅可比迭代次数已经超过最大迭代次数max1')enda,X=X;jX=X1',习题4.22.（1）取最大迭代次数Max1=100在MATLAB工作窗口输入程序>>A=[23 -1 -2;-1 10 -2;-1 -1 5];b=[1.7;8.3;4.2];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,100)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：a =-21 -8 -1jX =0.2159 1.0711 1.0974X =0.2159 1.0710 1.0973取最大迭代次数Max=5在MATLAB工作窗口输入程序>> A=[23 -1 -2;-1 10 -2;-1 -1 5];b=[1.7;8.3;4.2];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛请注意：雅可比迭代次数已经超过最大迭代次数max1a =-21 -8 -1jX =0.2159 1.0711 1.0974X =0.2152 1.0697 1.0959（2）取最大迭代次数Max1=100在MATLAB工作窗口输入程序>>A=[15 -1 -2;-1 10 -2;-1 -1 0.5];b=[7.2;8.3;4.2];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,100) 运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛取最大迭代次数Max1=5在MATLAB工作窗口输入程序>> A=[15 -1 -2;-1 10 -2;-1 -1 0.5];b=[7.2;8.3;4.2];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,5) 运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛（3）取最大迭代次数Max1=100在MATLAB工作窗口输入程序>> A=[10 -1 -2;-1 10 -2;-1 -1 5];b=[7.2;8.3;4.2];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,100) 运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：a =-8 -8 -1jX =1.1000 1.2000 1.3000X =1.0998 1.1998 1.2998取最大迭代次数Max1=5在MATLAB工作窗口输入程序>> A=[10 -1 -2;-1 10 -2;-1 -1 5];b=[7.2;8.3;4.2];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛请注意：雅可比迭代次数已经超过最大迭代次数max1a =-8 -8 -1jX =1.1000 1.2000 1.3000X =1.0951 1.1951 1.2941（4）取最大迭代次数Max1=100在MATLAB工作窗口输入程序>>A=[1 2 3;2 5 2;3 1 5];b=[14;18;20];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,100)运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛取最大迭代次数Max1=5在MATLAB工作窗口输入程序>> A=[1 2 3;2 5 2;3 1 5];b=[14;18;20];X0=[0 0 0]';X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛用高斯-塞德尔迭代定义解线性方程组的MATLAB主程序function X=gsdddy(A,b,X0,P,wucha,max1)D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1); dD=det(D);if dD==0disp('请注意：因为对角矩阵D奇异，所以此方程组无解.')elsedisp('请注意：因为对角矩阵D非奇异，所以此方程组有解.')iD=inv(D-L); B2=iD*U;f2=iD*b;jX=A\b;X=X0; [n m]=size(A);for k=1:max1X1= B2*X+f2; djwcX=norm(X1-X,P);xdwcX=djwcX/(norm(X,P)+eps);if (djwcX<wucha)|(xdwcX<wucha)breakelsek;X1';k=k+1;X=X1;endendif (djwcX<wucha)|(xdwcX<wucha)disp('请注意：高斯-塞德尔迭代收敛,此A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下： ')elsedisp('请注意：高斯-塞德尔迭代的结果没有达到给定的精度，并且迭代次数已经超过最大迭代次数max1,方程组的精确解jX和迭代向量X如下： ')X=X';jX=jX'endendX=X';D,U,L,jX=jX'习题4.33.（1）在MATLAB工作窗口输入程序>> A=[11 -1 -2;-1 10 -2;-1 -1 0.5];b=[7.2 ;8.3;4.2];X0=[0 0 0]';X=gsdddy(A,b,X0,inf,0.00001,100)运行后输出结果请注意：因为对角矩阵D非奇异，所以此方程组有解.请注意：高斯-塞德尔迭代收敛,此A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下：D =11.0000 0 00 10.0000 00 0 0.5000U =0 1 20 0 20 0 0L =0 0 01 0 01 1 0jX =15.8529 17.3941 74.8941X =15.8518 17.3928 74.8892（2）在MATLAB工作窗口输入程序>> A=[4 4 -5 7;2 -8 3 -2;4 5 -13 16;7 -2 2 3];b=[5;2;-1;21];X0=[0 0 0 0]';X=gsdddy(A,b,X0,inf,0.00001,100) 运行后输出结果请注意：因为对角矩阵D非奇异，所以此方程组有解.请注意：高斯-塞德尔迭代的结果没有达到给定的精度，并且迭代次数已经超过最大迭代次数max1,方程组的精确解jX和迭代向量X如下：jX =0.3446 0.7555 4.7894 3.5066D =4 0 0 00 -8 0 00 0 -13 00 0 0 3U =0 -4 5 -70 0 -3 20 0 0 -160 0 0 0L =0 0 0 0-2 0 0 0-4 -5 0 0-7 2 -2 0jX =0.34460.75554.78943.5066X =1.0e+26 *-0.7417-0.33740.07681.4545用谱半径判别超松弛迭代法产生的迭代序列的敛散性的MATLAB主程序function H=ddpbj(A,om)D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1); iD=inv(D-om*L);B2=iD*(om*U+(1-om)*D);H=eig(B2);mH=norm(H,inf);if mH>=1disp('请注意：因为谱半径不小于1，所以超松弛迭代序列发散,谱半径mH和B的所有的特征值H如下：')elsedisp('请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B的所有的特征值H如下：')endmH习题4.41.（1）当取ω=1.15时，在MATLAB工作窗口输入程序>> A=[7 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7];H=ddpbj(A,1.15)运行后输出结果请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B的所有的特征值H如下：mH =0.1608H =0.0715 + 0.1440i0.0715 - 0.1440i-0.1308 + 0.0498i-0.1308 - 0.0498i（2）当取ω=1.15时，在MATLAB工作窗口输入程序>> A=[7 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7];H=ddpbj(A,5)运行后输出结果请注意：因为谱半径不小于1，所以超松弛迭代序列发散,谱半径mH和B的所有的特征值H如下：mH =13.1892H =-13.1892 + 0.0000i-2.6969 + 0.0000i0.2460 + 2.6714i0.2460 - 2.6714i。