数值分析 上机实验
非线性方程求根
一、实验目的
1． 了解一般非线性方程的求根是比较复杂的事情：要讨论（或知道）它有无
实根，有多少实根；知道求近似根常用的几种方法，每种方法的特点是什么。
2． 用通过二分法（区间半分法）、不动点(也Picard)迭代法及Newton迭代（切线）法求其它非线性方程的根，并尽可能估计误差。
二、实验原理
三、实验程序
四、实验内容(对以下各题略加改动，但不要用这些题)
五、解答（按如下顺序提交电子版）
1.(程序)
二分法：
function x = erfen(a,b,TOL2) % 运用二分法求f(x)=0 % 求出的结果为 p % TOL2为距离误差 items = ceil(abs(log10((b-a)/TOL2)/log10(2))); syms t; f = t^3 - t - 1; if a>b disp('Method Failed !!!'); else count = 0; p = (a + b)/2; while (b - a) > TOL2 if subs(f,t,p)*subs(f,t,b) < 0 a = p; else b = p; end p = (a + b)/2; count = count + 1; if count > items break; end p = (a + b)/2; end x = p; End
不动点迭代法：
function [xc,num,eps] = budongdian(x0,phi,step) if nargin<2 phi = 10^(-6); end if nargin<3 step = 100; end preNum = x0; num = 0; eps =1; while eps > phi afterNum = g(preNum); eps = abs(g(afterNum) - g(preNum)); preNum = afterNum; num = num + 1; if num > step disp('超过迭代次数，可能不收敛！！！！'); break; end end xc = afterNum;
function y = g(t)
sym t;
y=(t+1)^(1/3);
NEWTON法：
function [x] = NEWTON(x0,m) TOL = 10^(-10); syms t; f = (t-1)^2; p0=x0; p = p0 - subs(diff(f),t,p0); while abs(p-p0) >= TOL for k = 1:m p0 = p; p = p0 - subs((f),t,p0)/subs(diff(f),t,p0); end end x = p0;
2.(运算结果)
二分法：题目：
>> x = erfen(0,5,10^(-6))
x =
1.3247
不动点迭代法：题目 ：
>> [xc,num,eps] = budongdian(2)
xc =
1.3247
num =
8
eps =
8.3221e-007
牛顿迭代法：题目：
>> [x] = NEWTON(100,10)
x =
1.0000
3.(拓展（方法改进、体会等）)
(1)使用二分法求解方程时，迭代次数items是由区间长度及精确度确定的。可以通过 ceil(abs(log10((b-a)/TOL2)/log10(2)))把最大迭代次数计算出来赋给迭代次数，使结果进行多次循环以求得精确解。
（2）不动点求根，首先要确定根的区间，在根附近选取已初值，经过多次迭代后，取满足误差判定的解作为最优解。
‚其中输入的函数应为迭代函数，且必须为收敛的（不动点迭代法的迭代函数仅为局部收敛，故不可任意给定初始值），且迭代函数应有连续的性质。
（3）牛顿迭代法求根，该函数必须连续可导，通过多次迭代可求出最优解。
矩阵的LU分解及在解线性方程组中的应用
一、实验目的
1． 借助矩阵理论进一步对消去法作分析，建立高斯消去法与矩阵因式分解的关系。
2． 会矩阵的紧凑格式的LU分解法、对称阵的分解法。
3． 会直接三角分解法线性方程组；会选列主元三角分解法解线性方程组。
二、实验原理
Gauss-Jordan消元法;初等变换。
三、实验程序
4. 将第2个程序改为选列主元三角分解法解方程组的程序。
四、实验内容
1． 求一个4阶矩阵的LU分解。
2． 用直接三角分解法，求一个4元线性方程组的解。
五、解答（按如下顺序提交电子版）
1.(程序)
(1)LU分解
function [L,A] = LU2(A) %Doolottle分解，A为输入的矩阵，L,A分别为下三角矩阵和上三角矩阵 [row_A,col_A] = size(A); L = zeros(row_A,row_A); for i = 1:row_A L(i,i) = 1; end n = size(A,1); for k = 1:n-1 for i = k+1:n for j = k + 1:n L(i,k) = A(i,k)/A(k,k); A(i,j) = A(i,j) - L(i,k)*A(k,j); end end end for i=1:n for j =1:n if(i>j) A(i,j) = 0; end end end
(2)直接三角解法
function [x]=zjT(A,b) %A代表Ax=B中的A； %b代表b； %x为最终求得的解。 [row_A,col_A] = size(A); B = zeros(row_A,col_A+1); for i = 1:row_A for j = 1:col_A B(i,j) = A(i,j); end end for j = 1:length(b) B(j,col_A+1) = b(j); end n = row_A; for k = 1:n-1 for i = k+1:n B(i,k) = B(i,k)/B(k,k); for j = k+1:n+1 B(i,j)= B(i,j)-B(i,k)*B(k,j); end end end for k = n:1:2 for i = k-1:1 B(i,k) = B(i,k)/B(k,k); B(i,n+1) = B(i,n+1) - B(i,k)*B(k,n+1); B(k,n+1) = B(k,n+1)/B(k,k); end end %求解 L 和 U L = zeros(n,n); U = zeros(n,n); for i = 1:n for j=1:n if i<=j U(i,j) = B(i,j); else L(i,j) = B(i,j); end if i == j L(i,j) = 1; end end end for k=1:n b(k) = B(k,col_A+1); end %回代计算 x = zeros(3,1); x(n,1) = b(n)./U(n,n); for i = 1:n-1 i=n-i; sum = 0; for k = i+1:n sum = sum + U(i,k)*x(k); end x(i,1) = (b(i) - sum)/U(i,i); end
2.(运算结果)
(1)LU分解
s =
2 3 4 5
4 8 13 16
6 17 36 44
8 22 65 86
>> [L,U]=LU2(s)
L =
1 0 0 0
2 1 0 0
3 4 1 0
4 5 6 1
U =
2 3 4 5
0 2 5 6
0 0 4 5
0 0 0 6
(2)直接三角分解
>> [x]=zjT(A,b)
x =
1
2
3
3.(拓展（方法改进、体会等）)
更优良的杜丽特尔解法：
function A=lu_doolittle(A) n=length(A); A(2:n,1)=A(2:n,1)/A(1,1); for k=2:n A(k,k:n)= A(k,k:n)-A(k,1:k-1)*A(1:k-1,k:n); A(k+1:n,k)=(A(k+1:n,k)-A(k+1:n,1:k-1)*A(1:k-1,k))/A(k,k); end
矩阵分解可以提高求解方程组的效率，还可以减小储存数据所需要的储存空间，这就是矩阵分解的优点。另外，如果主对角线元素较小，必须进行列主元交换，用绝对值最大的元素作为列主元。
线性方程组的迭代求解法
一、实验目的
1． 借助矩阵按模最大特征值，判断解方程组的Jacobi迭代法所得迭代序列的敛散性。
2． 会在Jacobi迭代法所得迭代序列收敛时，用修改后的Gauss-Seidel迭代法。
3． 会逐次超松驰迭代法。
二、实验原理
三、实验程序
四、实验内容
用上面前二种方法求解4元线性方程组的近似解，所选方程组尽可能可以用多种方法求得收敛解。
注：要注意判断迭代法收敛性，方法之一就是用程序求矩阵的按模最大特征值。
五、解答（按如下顺序提交电子版）
1.(程序)
（1）雅克比迭代法
function x=JCOB(A,b,x0) D = diag(diag(A)); B = D - A; k=0; while k<100 x=(inv(D)*B)*x0 + inv(D)*b; if norm(x-x0) < 0.0001 break; end x0=x; k=k+1; End
（2）高斯—赛德尔迭代法
function X=GS(A,b,X,m) %A,b为方程组的系数矩阵的伴随矩阵 %X为初值 %m为最大迭代次数 n=rank(A); for k=1:m for i=1:n sum=0; for j=i+1:n sum=sum+A(i,j)*X(j); end if i==1 X(i)=(b(i)-sum)/A(i,i); else sum1=0; for j=1:i-1 sum1=sum1+A(i,j)*X(j); end X(i)=(b(i)-sum1-sum)/A(i,i); end end end
2.(运算结果)
（1）
>> A=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4];b=[1;1;1;1];X0=[0;0;0;0];
>> JCOB(A,b,X0)
ans =
-0.9999
-0.9999
-0.9999
-0.9999
（2）
>> A=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4];b=[1;1;1;1];X0=[0;0;0;0];
>> GS(A,b,X0,100)
ans =
-1.0000
-1.0000
-1.0000
-1.0000
3.(拓展（方法改进、体会等）)
（3）超松弛迭代法
function [n,x]=sor22(A,b,X,nm,w,ww) %用超松弛迭代法求解方程组Ax=b %输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子 %输出：x为求得的方程组的解构成的列向量，n为迭代次数 n=1; m=length(A); D=diag(diag(A)); %令A=D-L-U,计算矩阵D L=tril(-A)+D; %令A=D-L-U,计算矩阵L U=triu(-A)+D; %令A=D-L-U,计算矩阵U M=inv(D-ww*L)*((1-ww)*D+ww*U); %计算迭代矩阵 g=ww*inv(D-ww*L)*b; %计算迭代格式中的常数项 %下面是迭代过程 while n<=nm x=M*X+g; %用迭代格式进行迭代 if norm(x-X,'inf')<w disp('迭代次数为');n disp('方程组的解为');x return; %上面：达到精度要求就结束程序，输出迭代次数和方程组的解 end X=x;n=n+1; end %下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果,(可能不是方程的解) disp('在最大迭代次数内不收敛!'); disp('最大迭代次数后的结果为');x
>> A=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4];b=[1;1;1;1];X0=[0;0;0;0];
>> [n,x]=sor22(A,b,X0,1000,0.0000001,1.3)
迭代次数为
n =
16
方程组的解为
x =
-1.0000
-1.0000
-1.0000
-1.0000
n =
16
x =
-1.0000
-1.0000
-1.0000
-1.0000
插值法
求上按5个等距节点确定的Lagrange, Newton插值多项式.
五、解答（按如下顺序提交电子版）
1.(程序)
（1）Lagrange插值多项式
function yy=lagrange1(x1,y1,xx) %本程序为Lagrange1插值，其中x1,y1 %为插值节点和节点上的函数值，输出为插值点xx的函数值， %xx可以是向量。 syms x; n=length(x1); for i=1:n t=x1; t(i)=[]; L(i)=prod((x-t)./(x1(i)-t)); % L向量用来存放插值基函数 end u=sum(L.*y1); p=simplify(u) % p是简化后的Lagrange插值函数（字符串） yy=subs(p,x,xx);
（2）Newton插值多项式
%Newton基本插值公式 %x为向量，全部的插值节点 %y为向量，差值节点处的函数值 function [c,d]=newton(x,y) n=length(x); d=zeros(n,n); d(:,1)=y'; format long; for j=2:n for k=j:n d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1)); end end c=d(n,n); for k=(n-1):-1:1 c=conv(c,poly(x(k))); m=length(c); c(m)=c(m)+d(k,k); end
2.(运算结果)
（1）拉格朗日法
>> yy=lagrange1(x,y,2.4)
p =
x^4
yy =
33.1776
（2）牛顿法
>> newton(x,y)
ans =
1.000000000000004 -0.000000000000020 0.000000000000031 -0.000000000000019 0.000000000000004
3.(拓展（方法改进、体会等）)
1.插值法是用几个已知的点去描述原函数的方法。
2.拉格朗日插值是通过建立拉格朗日基函数，从而描述原函数。拉格朗日基函数以及拉格朗日法具有代表性。它代表了一种解题的方法：设置特殊的基函数，通过求解基函数从而求解原方程。
最小二乘法
设有如下数据：
| x |
-3 |
-2 |
-1 |
0 |
1 |
2 |
3 |
| y |
-1.76 |
0.42 |
1.2 |
1.34 |
1.43 |
2.25 |
4.38 |
用3次多项式拟合这组数据.
解答（按如下顺序提交电子版）
1.(程序)
最小二乘法
function [p] = mincheng(x,y,n) %n代表拟合几次多项式 %x，y分别代表输入的数列 imax = max(x); imin = min(x); p=polyfit(x,y,n); x1=imin:0.05:imax; y1=polyval(p,x1); plot(x,y,'*r',x1,y1,'-b')
2.(运算结果)
>> x=[ -3 -2 -1 0 1 2 3]
x =
-3 -2 -1 0 1 2 3
>> y=[-1.76 0.42 1.2 1.34 1.43 2.25 4.38]
y =
-1.7600 0.4200 1.2000 1.3400 1.4300 2.2500 4.3800
>> mincheng(x,y,3)
ans =
0.1133 -0.0018 0.0035 1.3300
3.(拓展（方法改进、体会等）)
（1）最小二乘法实质上是在全数域上，求解一个使目标函数达到最小值的最优解。也就是让梯度向量等于零。（运筹学无约束非线性规划）
（2）最小二乘法也是最佳平方逼近法在离散的点上的应用，其法方程两端均为内积。
数值积分法
一、实验目的
许多工程技术和数学研究中要用到定积分，如果无法直接算不出精确值（如含在积分方程中的积分）或计算困难但可用近似值近似时，就用数值积分法方法加以解决。常用的算法有：复化梯形、辛甫生（Simpson）、柯特斯（Cotes）求积法; 龙贝格(Romberg)算法；高斯（Gauss）算法。
二、实验原理
三、实验程序
下面给出复化Ｓimpson求积法程序（梯形及柯特斯复化求积分程序可比照编制）：
四、实验内容
选一可精确算值的定积分，用复化的梯形法及复化Ｓimpson求积法作近似计算，并比较结果。
五、解答（按如下顺序提交电子版）
1.(程序)
(1)复化梯形公式
function S = Tixing(a,b,m) %复化辛普森公式求积公式程序 %S为最终的求得的近似值 %a,b为端点 %m为几等分 syms x; y = x/(4+x^2); T=0; h=(b-a)/m; for i=1:m-1 s = i*h + a; T=T+subs(y,x,s); end T=2*T+subs(y,x,a)+subs(y,x,b); S=T*h/2; S = vpa(S);
(2)复化Ｓimpson求积法
function S = simpson(a,b,m) %复化辛普森公式求积公式程序 %S为最终的求得的近似值 %a,b为端点 %m为几等分 syms x; y = x/(4+x^2); h = (b-a)/(2*m); s0 = subs(y,x,a)+subs(y,x,b); s1=0; s2=0; for i = 1:(2*m-1) l = a + i*h; if mod(i,2)==0 s2 = s2 + subs(y,x,l); else s1 = s1 + subs(y,x,l); end end s1 = h*(s0+4*s1+2*s2)/3; S = vpa(s1);
2.(运算结果)
(1)复化梯形公式
>> Tixing(0,1,8)
ans =
0.11140235452954801006222544000245
(2)复化Ｓimpson求积法
>> simpson(0,1,8)
ans =
0.1115718132526306492148607762825
3.(拓展（方法改进、体会等）)
（1）、数值积分的精度实质上是将估计公式中取为，如果满足求积公式，则说明积分公式至少具有m次精度；若求积公式在处不精确成立，则说明求积公式只具有m次精度。
（2）、数值积分实际上就是对一些无法求出原函数或者原函数比较复杂的积分进行计算的方法。这种方法来源于插值法，将目标函数在不同的节点进行插值，之后对插值多项式进行积分，从而近似目标值。
（3）、外推技巧：将两个精度较低的近似值通过线性组合的方式提高精度。我觉得这种方法非常精巧。将两个精度一致的近似值，如：，通过简单的线性组合来提高精度。其内涵是通过这种线性组合消掉了余项中的，提高精度
常微分方程的数值解法
一、实验目的
科学技术中常常要求解常微分方程的定解问题，所谓数值解法就是求未知函数在一系列离散点处的近似值。
二、实验内容
选一可求解的常微分方程的定解问题，分别用以上1, 4两种方法求出未知函数在
节点处的近似值，并对所求结果与分析解的（数值或图形）结果进行比较。
三、解答（按如下顺序提交电子版）
1.(程序)
（1）欧拉法
function y = oula(a,b,n,y0) %欧拉法 %a,b位区间端点 %n为区间分的分数 %y0代表初值 h = (b-a)/n; t = a; y = y0; for i = 1:n y = y + h*f(t,y); t = a + i*h; end function p = f(u,v) p = v - (2*u)/v;
（2）欧拉两步法
%两步欧拉法 function y = danbufa(a,b,n,y0,y1) %a,b位区间端点 %n为区间分的分数 %y0,y1代表初值 h = (b-a)/n; t = a; y(1) = y0; y(2) = y1; for i = 2:n-1 y(i+1) = y(i-1) + h*f(t,y(i)); t = a + i*h; end function p = f(u,v) p = v - (2*u)/v;
2.(运算结果)
（1）欧拉法
>> y = oula(0,1,10,1)
y =
1.7848
（2）两部欧拉法
>> y = danbufa(0,1,10,1,1.1)
y =
1.0000 1.1000 1.1100 1.1750 1.1764 1.2246 1.2172 1.2477 1.2298 1.2406
总结： 这学期学习了数值分析的课程，学会了很多算法，虽然以后可能不会从事数学工作，但还是好好学了一遍，本门课程真的很有用！
浙公网安备 33010602011771号