《计算方法》实验报告指导教师:学院:班级:团队成员:一、题目例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解k x ,并使其满足8110k k x x ---<原理:在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ,过曲线()y f x =的点()()00x f x ,引切线()()()1000:'l y f x f x x x =+-其与x 轴相交于点:()()0100 'f x x x f x =-,进一步,过曲线()y f x =的点()()11x f x , 引切线()()()2111: 'l y f x f x x x =+-其与x 轴相交于点:()()1211 'f x x x f x =-如此循环往复,可得一列逼近方程()0f x =精确解*x 的点01k x x x ,,,,,其一般表达式为:()()111 'k k k k f x x x f x ---=-该公式所表述的求解方法称为Newton 迭代法或切线法。
程序:function y=f(x)%定义原函数y=x^3-x-1;endfunction y1=f1(x0)%求导函数在x0点的值syms x;t=diff(f(x),x);y1=subs(t,x,x0);endfunction newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数while abs(x1-x0)>=tolx0=x1;x1=x0-f(x0)/f1(x0);k=k+1;endfprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1);fprintf('迭代次数为k=%d\n',k);end结果:二、题目例3.7试利用Jacobi 迭代公式求解方程组1234451111101112115181111034x x x x ----⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪---⎝⎭⎝⎭⎝⎭ 要求数值解()k X 满足()4210k --≤X X ,其中T (1,2,3,4)=X 为方程组的精确解。
原理:将线性方程组的系数矩阵111212122212n n n n nn a a a a a a a a a ⎛⎫ ⎪⎪= ⎪ ⎪⎝⎭A 分解为=++A L D U ,其中1122=(,,,)nn diag a a a D ,3132,12112000000000=0n n n n a a a a a a -⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭L ,121312321,00=00000000n n n n a a a a a a -⎛⎫ ⎪⎪⎪ ⎪ ⎪ ⎪⎝⎭U 当对角阵D 可逆时,方程组=AX b 可等价地写成()11 --=-++X D L U X D b ,据此可得Jacobi 迭代公式为()()()111 0,1,k kk +--=-++=X D L U X D b, ,其中()()()()()T12=,k k k n nk x x x ∈XR ,,,该迭代公式也可以写成如下的分量形式1(()11),,1,2,,1nk k i ij j j j i ii ixb a i n a x ++=≠⎛⎫= ⎪⎭=⎝-∑程序:function jacobi()%输入矩阵A 、b 、精度tol%A=[5 -1 -1 -1; -1 10 -1 -1; -1 -1 5 -1; -1 -1 -1 10]; b=[-4 12 8 34]'; tol=10^-4;% A=input('系数矩阵A=');%输入矩阵A % b=input('矩阵b= ');%输入矩阵b% tol=input('精度要求tol=');%输入精度tol X=inv(A)*b; [n,~]=size(A); D=diag(diag(A)); L=tril(A)-D; U=triu(A)-D; X0=zeros(n,1);X1=inv(D)*b-inv(D)*(L+U)*X0; X0=X1; X2=X-X0; k=1;while abs(norm(X2,2))>=tolX1=inv(D)*b-inv(D)*(L+U)*X0; X0=X1; X3=X-X0; %收敛性判断%if norm(X3,2)>norm(X2,2) break; elseX2=X3; k=k+1; end end% 结果输出%if X2~=X3fprintf('应用jacobi迭代法不收敛\n');elsefprintf('\n应用jacobi迭代公式迭代 k=%d次后\n可得满足精度要求的数值解 X(%d)=\n',k,k);for i=1:nfprintf(' %1.15f\n',X1(i));endfprintf('且其满足计算精度norm(abs(X-X(%d)),2)=%1.15g<%.1e\n',k,norm(abs(X-X1),2),tol);end结果:三、题目例3.8试利用Gauss Seidel -迭代公式求解方程组1234451111101112115181111034x x x x ----⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪---⎝⎭⎝⎭⎝⎭ 要求数值解()k X 满足()4210k --≤X X ,其中T (1,2,3,4)=X 为方程组的精确解。
原理:将线性方程组的系数矩阵111212122212n n n n nn a a a a a a a a a ⎛⎫ ⎪⎪= ⎪ ⎪⎝⎭A 分解为=++A L D U ,其中1122=(,,,)nn diag a a a D ,3132,12112000000000=0n n n n a a a a a a -⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭L ,121312321,00=00000000n n n n a a a a a a -⎛⎫ ⎪⎪⎪ ⎪ ⎪ ⎪⎝⎭U 当对角阵D 可逆时,方程组=AX b 可等价地写成=--+DX LX UX b ,由此可构造迭代格式,(1)(1)()k k k b ++=--+DXLX UX据此可得Gauss Seidel -迭代公式为()111()(1)0,1,k k k -++-=-+++=X D L UX D L b (), ,该公式也可以写成如下的分量形式11())11(()11),1,2(,,i nk k k i i ij j ij j j j i ii x b a x a x i n a -++==+=--=∑∑程序:function Gauss_Seidel()%输入矩阵A 、b 、精度tol%A=[5 -1 -1 -1; -1 10 -1 -1; -1 -1 5 -1; -1 -1 -1 10]; b=[-4 12 8 34]'; tol=10^-4;% A=input('系数矩阵A=');%输入矩阵A % b=input('矩阵b=');%输入矩阵b% tol=input('精度要求tol=');%输入精度tol X=inv(A)*b; [n,~]=size(A); D=diag(diag(A)); L=tril(A)-D; U=triu(A)-D; X0=zeros(n,1);X1=inv(D+L)*b-inv(D+L)*U*X0; X0=X1; X2=X-X0; k=1;while abs(norm(X2,2))>=tolX1=inv(D+L)*b-inv(D+L)*U*X0; X0=X1; X3=X-X0; %收敛性判断%if norm(X3,2)>norm(X2,2) break; elseX2=X3; k=k+1; end end% 结果输出%if X2~=X3fprintf('应用Gauss_Seidel迭代法不收敛\n');elsefprintf('\n应用Gauss_Seidel迭代公式迭代%d次后\n可得满足精度要求的数值解 X(%d)=\n',k,k);for i=1:nfprintf(' %1.15f\n',X1(i));endfprintf('且其满足计算精度norm(abs(X-X(%d)),2) =%1.15g<%.1e\n',k,norm(abs(X-X1),2),tol);end结果:四、题目例4.1给定函数()(1cos )f x x x =+及插值节点原理:设函数()f x 在区间[],a b 上有定义,其在1n +个互异点[],i x a b ∈ 2,,n .代数插值问题,所求得的(n P x 0n > ,1,,0,1,2,,0,i ji j n i j=⎧=⎨≠⎩,.程序:function y=f(x)%定义原函数y=x*(1+cos(x));endfunction y1=maxM(a,b)%输入区间[a,b]的上下限%syms xy=f(x);y5=diff(y,5);%求n+1阶导y1=max(abs(subs(y5,x,[a:0.01:b])));%求导函数的绝对值在[a,b]区间上的最大值endfunction lagrange_interpolation(x)%输入待求点x值矩阵x% 给定插值节点x0[N]%x0=[0,pi/8,pi/4,3*pi/8,pi/2];% x0=input('x0=');%输入插值节点矩阵x0n=length(x0);m=length(x);% 计算插值节点所对应的函数值%for i=1:ny0(i)=f(x0(i));end% lagrange插值法%for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endw=1.0;% 结果输出%for i=1:mfor j=1:nw=w*(x(i)-x0(j));endfprintf('f(x(%d))的近似值为 %1.15g\n',i,y(i));fprintf(' 误差估计<=%1.15e\n',maxM(x0(1),x0(n))/prod(1:n)*abs(w));%调用MaxM函数 fprintf(' 绝对误差为 %1.15e\n',abs(f(x(i))-y(i))); endend结果:五、题目例4.3已知()f x =012340.35,0.50,0.65,0.80,0.95,x x x x x =====构造4次Newton 插值多项式,由此计算(0.7)f 的逼近值,并指出其绝对误差。