《数值分析matlab程序.pdf》由会员分享，可在线阅读，更多相关《数值分析matlab程序.pdf（46页珍藏版）》请在三一文库上搜索。
1、数值分析（ matlab程序） 曹德欣 曹璎珞 第一章绪论 数值稳定性程序，计算P11 试验题一积分 function try_stable global n global a a=input(a=); N = 20; I0 = log(1+a)-log(a); I = zeros(N,1); I(1) = -a*I0+1; for k = 2:N I(k) = - a*I(k-1)+1/k; end II = zeros(N,1); if a=N/(N+1) II(N) = (1+2*a)/(2*a*(a+1)*(N+1); else II(N) =(1/(a+1)*(N+1)+1/N)/2
2、; end for k = N:-1:2 II(k-1) = ( - II(k) +1/k) / a; end III = zeros(N,1); for k = 1:N n = k; III(k) = quadl(f,0,1); end clc fprintf(n 算法 1 结果算法 2 结果精确值 ) for k = 1:N, fprintf(nI(%2.0f) %17.7f %17.7f %17.7f,k,I(k),II(k),III(k) end function y = f(x) global n global a y = x.n./(a+x); return 第二章非线性方程求解
3、下面均以方程y=x4+2*x2-x-3为例： 1、二分法 function y=erfen(a,b,esp) format long if narginesp if fun(a)*fun(c)=1.0e-4) % 对增广矩阵实施行交换 end if Aug(k,k)=0, error(这是奇异矩阵!) % 程序遇到error 会中断执行并显示其中的提示内 容 end % 把增广矩阵消元成为上三角 for p = k+1:n mult = Aug(p,k)/Aug(k,k); % 消元乘子 Aug(p,k:n+1) = Aug(p,k:n+1) - mult * Aug(k,k:n+1); en
4、d end % 解上三角方程组 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k) = ( b(k) - A(k,k+1:n)*x(k+1:n) ) / A(k,k); end 3.4 % exp3_4.m - 病态矩阵试验 % Hilbert 矩阵是严重的病态的矩阵(见 P56),求解病态方程组是相当困难 % 一般的方法都会失效,例如列主元Gauss 消去法 ,要采取特殊的方法和技术 % 例如(1) pcg (预处理共轭梯度法,只适用正定方程组), % (2) gmres (广义最残差法,适用一般
5、的大型疏矩阵) n = 15; A = hilb(n); % n 阶 Hilbert 矩阵 x = ones(n,1); % 指定精确解全是1 b = A*x; c = cond(A,inf); % 求条件数 ,范数取最大值范数 x1 = Ab; % 用 Gauss列主元消去法求解或者调用你编的x1 = gauss(A,b) x2 = pcg(A,b); % 用 pcg 法求解 (A 必须是正定矩阵) x3=gmres(A,b); % gmres 法求解结果 % 显示结果 clc fprintf(A 的条件数= %fn,c) fprintf(nGauss 法求解结果pcg 法求解结果gmres
6、 法求解结果 ) for i = 1:n fprintf(n %6.2f %6.2f %6.2f,x1(i),x2(i),x3(i) end % 注：范德蒙 (Vandermonde)矩阵也是严重病态矩阵 % A = vander(V), 其中 V 是一向量 % 对于维数较大的V,看一下范德蒙矩阵的条件数 3.5 % exp3_5.m - SOR 迭代法收敛速度受松驰因子的影响试验 function try_sor_and_relaxation % 参见P64 例 8 和 P65 表 3-4 A = 2 -1 0 0 -1 2 -1 0 0 -1 2 -1 0 0 -1 2; b = 1 0
7、1 0; tol = 1e-6; maxiter = 1000; P = 1 0 1 0; clc fprintf(* SOR迭代 :松驰因子与迭代次数的关系*n), fprintf(n 松驰因子迭代次数 ) fprintf(n-) for W = 1:0.1:1.9 X,iternum = sor(A,b,tol,maxiter,P,W); fprintf(n %3.2f %4.0f,W,iternum) end % 通过以上结果 ,你认为最佳松驰因子大致为多少? % - 说明- % 对于对称正定三对角矩阵(上面矩阵就是),其最佳松驰因子是可以求得的 % Wopt = 2 / ( 1+sqr
8、t(1-rho2) ) % 其中 rho 是 Jacobi 迭代矩阵Mj 的谱半径 % 下面求之 ,应该与我们实验的结果差不多吧,看看 D = diag(diag(A); Mj = eye(4) - inv(D)*A; rho = max(abs(eig(Mj); Wopt = 2 / ( 1+sqrt(1-rho2) ); fprintf(n-) fprintf(n 理论上最佳松驰因子是Wopt = %3.2f,Wopt) % - SOR 迭代法- function X,iternum = sor(A,B,tol,maxiter,P,W) % X,iternum = sor(A,B,tol,
9、maxiter,P,W) SOR迭代法解线性方程组AX=B % 输入- tol: 相对残向量的容差,当 norm(B-AX)/norm(B) epsilon % 找绝对值最大元素所在行与列(p,q) m1 p=max(abs(A-diag(diag(A); m2 q = max(m1); p=p(q); % 计算 Givens 矩阵 , 见 P197 (8-21) d = ( A(q,q)-A(p,p) ) / ( 2*A(p,q) ); if (d = 0) , t = 1/(d+sqrt(d*d+1); else t = 1/(d-sqrt(d*d+1); end c = 1/sqrt(1+t2); s = c*t; G = c,s;-s,c; % 消零 , 见 P197-198 (8-22)-(8-27) A(p q,:) = G*A(p q,:); A(:,p q) = A(:,p q)*G; Q(:,p q) = Q(:,p q)*G; alpha = norm(A-diag(diag(A),fro) / n; end D = diag(A); % -