贵州大学数值分析上机实验答案
数值分析上机作业(总)
数值分析上机实验一、解线性方程组直接法（教材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的数。
数值分析上机答案computer中科院
源程序如下：function[x_jac,x_gauss,x_sor,A,b]=diedai(n,w,e,N,Ain,bin) if nargin<3e=1e-5;N=1000;endif nargin<4N=1000;endif nargin<2error('²ÎÊý²»¹»');endA=hilb(n);b=A*ones(n,1);if nargin>4if nargin==5error('ÇëÊäÈëb');elseA=Ain;b=bin;endend%%Ñſ˱ȵü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ediedai_n=diedai_n+1;if diedai_n>Ndisp('Ñſ˱ȵü´úδ´ïµ½¾«¶È');break;endx(:,1)=x(:,2);for i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),1)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i);endendendx_jac.x=x(:,2);x_jac.n=diedai_n-1;x_jac.error=max(abs(b-A*x(:,2)));%%¸ß˹µü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('¸ß˹¡ªÈüµÂ¶ûµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i); endendendx_gauss.x=x(:,2);x_gauss.n=diedai_n-1;x_gauss.error=max(abs(b-A*x(:,2)));%%SORµü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('SORµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=x(i,1)+w*(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,i:n)* x(i:n,1))/A(i,i);elsex(i,2)=x(i,1)+w*(b(i)-A(i,i:n)*x(i:n,1))/A(i,i);endendendx_sor.x=x(:,2);x_sor.n=diedai_n-1;x_sor.error=max(abs(b-A*x(:,2)));end运行结果：由于SOR的松弛因子取1时与高斯-赛德尔迭代相同，故在矩阵n取8,10时省去了SOR迭代松弛因子取1的步骤。
数值分析实验题答案
数值分析实验报告姓名：院系：能源学院热能工程学号：2014年4月习题一实验3.2编制正交化多项式最小二乘拟合程序，并用于求解3次多项式最小二乘拟合为基的多项式最小问题，作拟合曲线的图形，计算平方误差，并与以函数{}0n k kx=二乘拟合的结果作比较表1x-1.0-0.50.00.5 1.0 1.5 2.0iy-4.447-0.4520.5510.0480.4470.549 4.552 i首先使用以函数{}0n k k为基的多项式最小二乘拟合，代码如下：x=然后使用正交化多项式方法作最小二乘拟合并画图，代码如下：拟合得到的图形如下（图1）：从图形来看，二者与数据点都很吻合。
计算结果为：delta1=2.1762e-05delta2=4.4701e-04为基的多项式拟合精度更高。
可以看出对于3次多项式以{}0n k kx=图1习题二 实验4.2分别用复化Simpson 公式与变步长Simpson 公式计算，要求绝对误差限为71=102ε-⨯，输出每种方法所需的节点数和积分近似值并分析比较结果。
（1）6220()10x x x dx -+⎰ （2）10⎰ （3）6220()10x x x dx -+⎰对于复化Simpson 公式，使用事前误差估计法得到所需计算节点数，有以下误差估计公式：4(4)()()()(),(,)1802n b a h R f f a b ηη-=-∈对（1）式有：(4)22max ()36144x f x x ===0.0266h ≤75.2b aN h -≥=取76N =，计算节点数为21153N N =+= 对（2）（3）式由于其4阶导数值分布极不均匀，用最大值来估计所需计算节点数造成很大浪费，尝试多次后分别取153N =和1091N =代码如下：计算结果如下：（1）计算结果用复化Simpson公式计算：节点数：153近似值：1.161904777用变步长Simpson公式计算：节点数：77近似值：1.161904766标准值：1.161904762（2）计算结果用复化Simpson公式计算：节点数：153近似值：0.400000049用变步长Simpson公式计算：节点数：33近似值：0.400000069标准值：0.400000000（3）计算结果用复化Simpson公式计算：节点数：1091近似值：23.812135331用变步长Simpson公式计算：节点数：145近似值：23.812135297标准值：23.812135292从以上计算结果可以看出，变步长simpson公式所需节点数明显减少，因为3个函数的4阶导数在积分区间内分布都部均匀，（2）（3）式更为严重，在先范围区间内导数值远大于其他区间，只需要在这些区间增加节点数就可以达到指定精度，而Simpson公式需要在全体积分限内采用较小间距才满足条件。
数值分析实验报告 课后答案_【khdaw_lxywyl】
包含工大硕士生数值分析课上机试验所有试题及答案 实验一：非线性方程求解 程序1：二分法： syms f x;f=input('请输入f(x)=');A=input('请输入根的估计范围[a,b]='); e=input('请输入根的误差限e='); while (A(2)-A(1))>e c=(A(1)+A(2))/2; x=A(1); f1=eval(f); x=c;f2=eval(f); if (f1*f2)>0 A(1)=c; elseA(2)=c; end endc=(A(1)+A(2))/2;fprintf('c=%.6f\na=%.6f\nb=%.6f\n',c,A)用二分法计算方程：1．请输入f(x)=sin(x)-x^2/2请输入根的估计范围[a,b]=[1,2] 请输入根的误差限e=0.5e-005 c=1.404413 a=1.404411 b=1.4044152．请输入f(x)=x^3-x-1请输入根的估计范围[a,b]=[1,1.5] 请输入根的误差限e=0.5e-005 c=1.324717 a=1.324715 b=1.324718程序2：newton 法： syms f x;f=input('请输入f(x)=');w ww .k hd aw .c o m课后答案网df=diff(f);x0=input('请输入迭代初值x0='); e1=input('请输入奇异判断e1='); e2=input('请输入根的误差限e2='); N=input('请输入迭代次数限N='); k=1;while (k<N) x=x0; if abs(eval(f))<e1fprintf('奇异!\nx=%.6f\n 迭代次数为:%d\n',x0,k) break elsex1=x0-eval(f)/eval(df); if abs(x1-x0)<e2fprintf('x=%.6f\n 迭代次数为:%d\n',x1,k) break elsex0=x1; k=k+1; end end end if k>=Nfprintf('失败\n') end用newton 法计算方程： 1．请输入f(x)=x*exp(x)-1请输入迭代初值x0=0.5请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=10 x=0.567143 迭代次数为:42．请输入f(x)=x^3-x-1请输入迭代初值x0=1请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=10w ww .k hd aw .c o m课后答案网x=1.324718 迭代次数为:53.1：请输入f(x)=(x-1)^2*(2*x-1)请输入迭代初值x0=0.45 请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=10 x=0.500000 迭代次数为:43.2：请输入f(x)=(x-1)^2*(2*x-1)请输入迭代初值x0=0.65 请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=10 x=0.500000 迭代次数为:93.3：请输入f(x)=(x-1)^2*(2*x-1)请输入迭代初值x0=0.55 请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=10 x=0.500000 迭代次数为:4程序3：改进的newton 法： syms f x;f=input('请输入f(x)='); df=diff(f);x0=input('请输入迭代初值x0='); e1=input('请输入奇异判断e1='); e2=input('请输入根的误差限e2='); N=input('请输入迭代次数限N='); k=1;while (k<N) x=x0; if abs(eval(f))<e1fprintf('奇异!\nx=%.6f\n 迭代次数为:%d\n',x0,k) break elsew ww .k hd aw .c o m课后答案网x1=x0-2*eval(f)/eval(df); if abs(x1-x0)<e2fprintf('x=%.6f\n 迭代次数为:%d\n',x1,k) break elsex0=x1; k=k+1; end end end if k>=Nfprintf('失败\n') end用改进的newton 法计算方程： 1．请输入f(x)=(x-1)^2*(2*x-1)请输入迭代初值x0=0.55 请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=10 失败2．请输入f(x)=(x-1)^2*(2*x-1)请输入迭代初值x0=0.55 请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=20 失败3．请输入f(x)=(x-1)^2*(2*x-1)请输入迭代初值x0=0.55 请输入奇异判断e1=0.1e-010 请输入根的误差限e2=0.5e-005 请输入迭代次数限N=100 失败实验二：Gauss 列主元消去法 程序1：Gauss 列主元消去法A=input('请输入线性方程组的增广矩阵A='); n=length(A)-1; x=zeros(n,1);w ww .k hd aw .c o m课后答案网aa=zeros(n,1); for j=1:n for i=1:(n+1)AA(j,i)=abs(A(j,i)); end endfor k=1:(n-1) for i=k:naa(i-(k-1))=AA(i,k); end for i=k:nif AA(i,k)==max(aa) break end end if AA(i,k)==0 breakfprintf('方程组系数矩阵奇异\n'); else for j=k:(n+1) jh=A(i,j); A(i,j)=A(k,j); A(k,j)=jh; end endfenzi=A(k,k); for j=k:(n+1)A(k,j)=A(k,j)/fenzi; end for p=(k+1):n jj=A(p,k); for j=k:(n+1)A(p,j)=A(p,j)-jj*A(k,j); end end endif k==(n-1)x(n)=A(n,(n+1))/A(n,n); for i=(n-1):(-1):1w ww .k hd aw .c o m课后答案网he=0; for j=(i+1):nhe=he+A(i,j)*x(j); endx(i)=A(i,(n+1))-he; end end x用Gauss 列主元消去法解方程组：1．请输入线性方程组的增广矩阵A=[1e-008,2,3,1;-1,3.172,4.623,2;-2,1.072,5.643,3]x =-0.4653 -0.0700 0.38002．请输入线性方程组的增广矩阵A=[4,-2,4,10;-2,17,10,3;-4,10,9,-7];x =2.9464 0.6071 -0.14293．请输入线性方程组的增广矩阵A=[0.3e-020,1,0.7;1,1,0.9]x =0.20000.7000程序2：不选主元的高斯消去法A=input('请输入线性方程组的增广矩阵A='); n=length(A)-1; x=zeros(n,1); for k=1:(n-1) if A(k,k)==0 breakfprintf('方程组不能用普通的高斯消去法解\n'); elsefenzi=A(k,k); for j=k:(n+1)A(k,j)=A(k,j)/fenzi; endw ww .k hd aw .c o m课后答案网for p=(k+1):n jj=A(p,k); for j=k:(n+1)A(p,j)=A(p,j)-jj*A(k,j); end endx(n)=A(n,(n+1))/A(n,n); for i=(n-1):(-1):1 he=0;for j=(i+1):nhe=he+A(i,j)*x(j); endx(i)=A(i,(n+1))-he; end end end x用不选主元的Gauss 消去法解方程组：1．请输入线性方程组的增广矩阵A=[4,-2,4,10;-2,17,10,3;-4,10,9,-7];x =2.9464 0.6071 -0.14292．请输入线性方程组的增广矩阵A=[1e-008,2,3,1;-1,3.172,4.623,2;-2,1.072,5.643,3];x =-0.4653 -0.07000.38003．请输入线性方程组的增广矩阵A=[0.3e-020,1,0.7;1,1,0.9]x =0 0.7000实验三:Runge 现象的产生和克服 程序1：lagrange 多项式插值 syms f x p dp lx L; f=1/(1+25*x^2);N=input('请输入插值节点数N=');w ww .k hd aw .c o m课后答案网xx=-1:2/N:1; p=1; L=0;ff=zeros(1,length(xx)); for i=1:(N+1) x=xx(i); ff(i)=eval(f); syms x;p=p*(x-xx(i)); enddp=diff(p); for j=1:(N+1) x=xx(j); k=eval(dp); syms x;lx=p/((x-xx(j))*k); L=L+lx*ff(j); endaa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96]; for i=1:length(aa) x=aa(i);S(i)=eval(L); fff(i)=eval(f); end e=0;for i=1:length(aa)e=e+(S(i)-fff(i))^2; ende=sqrt(e/(20*21));fprintf('插值偏差为e=%.6f\n',e) ezplot(f,[-1,1]) hold onezplot(L,[-1,1]) hold onplot(xx,ff,'*') hold on plot(aa,S,'o') hold off程序2：分段线性插值 syms f x p lx; f=1/(1+25*x^2);N=input('请输入插值节点数N='); xx=-1:2/N:1; p=1; L=0;ff=zeros(1,length(xx)); for i=1:(N+1) x=xx(i); ff(i)=eval(f); end syms xw ww .k hd aw .c o m课后答案网for i=1:N for j=1:(N+1) if j==ilx(i,j)=(x-xx(i+1))/(xx(i)-xx(i+1)); else if j==i+1lx(i,j)=(x-xx(i))/(xx(i+1)-xx(i)); elselx(i,j)=0; end end end end p=lx*ff';aa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96]; for i=1:length(aa) x=aa(i); for j=1:N+1 if x<xx(j) break end endS(i)=eval(p(j-1)); fff(i)=eval(f); end e=0;for i=1:length(aa)e=e+(S(i)-fff(i))^2; ende=sqrt(e/(20*21));fprintf('插值偏差为e=%.6f\n',e) ezplot(f,[-1,1]) hold onxxx=(-1:0.01:1); for i=1:length(xxx) x=xxx(i); for j=1:N+1 if x<xx(j) break end endSS(i)=eval(p(j-1)); endplot(xxx,SS,'r') hold onplot(xx,ff,'*') hold on plot(aa,S,'o') hold off程序3：三转角插值法w ww .k hd aw .c o m课后答案网syms f x df s s1 s2 s3 s4; f=1/(1+25*x^2); df=diff(f);N=input('请输入插值节点数N='); h=2/N;xx=-1:2/N:1; p=1; L=0;ff=zeros(1,length(xx)); for i=1:(N+1) x=xx(i); ff(i)=eval(f); dff(i)=eval(df); end syms x for i=1:Ns1=(x-xx(i+1))^2*(h+2*(x-xx(i)))*ff(i)/h^3; s2=(x-xx(i))^3*(h+2*(xx(i+1)-x))*ff(i+1)/h^3; s3=(x-xx(i+1))^2*(x-xx(i))*dff(i)/h^2; s4=(x-xx(i))^2*(x-xx(i+1))*dff(i+1)/h^2; s(i)=s1+s2+s3+s4; endaa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96]; for i=1:length(aa) x=aa(i); for j=1:N+1 if x<xx(j) break end endS(i)=eval(s(j-1)); fff(i)=eval(f); end e=0;for i=1:length(aa)e=e+(S(i)-fff(i))^2; ende=sqrt(e/(20*21));fprintf('插值偏差为e=%.6f\n',e) ezplot(f,[-1,1]) hold onxxx=(-1:0.01:1); for i=1:length(xxx) x=xxx(i); for j=1:N+1 if x<xx(j) break end endSS(i)=eval(s(j-1)); endw ww .k hd aw .c o m课后答案网plot(xxx,SS,'r') hold onplot(xx,ff,'*') hold on plot(aa,S,'o') hold off程序4：三弯矩插值法： syms f x ddf s; f=1/(1+25*x^2); ddf=diff(diff(f));N=input('请输入插值节点数N='); h=2/N;xx=-1:2/N:1; p=1; L=0;ff=zeros(1,length(xx)); for i=1:(N+1) x=xx(i); ff(i)=eval(f);ddff(i)=eval(ddf); end syms x for i=1:NA=(ff(i+1)-ff(i))/h-h*(ddff(i+1)-ddff(i))/6; B=ff(i)-h^2*ddff(i)/6;s(i)=(xx(i+1)-x)^3*ddff(i)/(6*h)+(x-xx(i))^3*ddff(i+1)/(6*h)+A*(x-xx(i))+B; endaa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96]; for i=1:length(aa) x=aa(i); for j=1:N+1 if x<xx(j) break end endS(i)=eval(s(j-1)); fff(i)=eval(f); end e=0;for i=1:length(aa)e=e+(S(i)-fff(i))^2; ende=sqrt(e/(20*21));fprintf('插值偏差为e=%.6f\n',e) ezplot(f,[-1,1]) hold onxxx=(-1:0.01:1); for i=1:length(xxx) x=xxx(i);w ww .k hd aw .c o m课后答案网for j=1:N+1 if x<xx(j) break end endSS(i)=eval(s(j-1)); endplot(xxx,SS,'r') hold onplot(xx,ff,'*') hold on plot(aa,S,'o') hold off图1：观察Runge 现象图2：分段线性插值w ww .k hd aw .c o m课后答案网图3：三转角插值：w ww .k hd aw .c o m课后答案网图4：三弯矩插值：w ww .k hd aw .c o m课后答案网实验四：多项式最小二乘法程序1：要求给出插值次数的多项式最小二乘法 syms x f;xx=input('请输入插值节点 as [x1,x2...]\n'); ff=input('请输入插值节点处对应的函数值 as [f1,f2...]\n'); m=input('请输入要求的插值次数m='); n=length(xx); for i=1:(m+1) syms fai x; fai=x^(i-1); for j=1:n x=xx(j);H(i,j)=eval(fai); end endA=ff*(H)'*inv(H*(H)'); syms x; f=0; for i=1:(m+1)w ww .k hd aw .c o m课后答案网f=f+A(i)*x^(i-1); endplot(xx,ff,'*') hold onezplot(f,[xx(1),xx(n)])用此程序作出课本上的两道题： 1. 请输入插值节点 as [x1,x2...] [0 0.9 1.9 3 3.9 5.0]请输入插值节点处对应的函数值 as [f1,f2...] [0 10 30 50 80 110]请输入要求的插值次数m=2 图：2.请输入插值节点 as [x1,x2...] [0 0.9 1.9 33.9 5.0]请输入插值节点处对应的函数值 as [f1,f2...] [19.0 32.3 49.0 73.3 97.8] 请输入要求的插值次数m=2 图：w ww .k hd aw .c o m课后答案网程序2：自适应多项式最小二乘法 syms x f nihef nihe;xx=input('请输入插值节点 as [x1,x2...]\n'); ff=input('请输入插值节点处对应的函数值 as [f1,f2...]\n'); n=length(xx); for m=1:(n-1) for i=1:(m+1) syms fai x; fai=x^(i-1); for j=1:nx=xx(j);H(i,j)=eval(fai); end endA=ff*(H)'*inv(H*(H)'); syms x f; f=0; for i=1:(m+1)f=f+A(i)*x^(i-1);w ww .k hd aw .c o m课后答案网endnihef(m)=f; for i=1:n x=xx(i);niheff(i)=eval(f); ee(i)=ff(i)-niheff(i); end e(m)=0; for i=1:(n-1)e(m)=(ee(i))^2+e(m); ende(m)=sqrt(e(m)/(n*(n-1))); eee(m)=1/e(m); endfor k=1:n if eee(k)==max(eee); break end end syms x;nihe=nihef(k); emcino=e(k); plot(xx,ff,'*') hold onezplot(nihe,[xx(1),xx(n)])fprintf('插值误差为e=%.6f',emcino) 用此程序作实验中要求的题目： 3.请输入插值节点 as [x1,x2...] [3 4 5 6 7 8 9]请输入插值节点处对应的函数值 as [f1,f2...] [2.01 2.98 3.50 5.02 5.47 6.02 7.05] 插值误差为e=0.000010 图：w ww .k hd aw .c o m课后答案网实验五：龙贝格积分法 程序：龙贝格积分： syms f x;f=input('请输入积分函数f='); A=input('请输入积分区间[a,b]='); e=input('请输入误差限e='); x=A(1); fa=eval(f); x=A(2); fb=eval(f);T(1)=(A(2)-A(1))*(fa+fb)/2; m=1;x=(A(1)+A(2));t=T(1)/2+(A(2)-A(1))*eval(f)/2; while abs(t-T(1))>e t=T(1); new=0; for i=1:(2^m)x=A(1)+(2*i-1)*(A(2)-A(1))/(2^(m+1));w ww .k hd aw .c o m课后答案网new=new+eval(f); endT(m+1)=T(m)/2+(A(2)-A(1))*new/(2^(m+1)); for i=m:(-1):1 L(m,i)=T(i); end for p=m:(-1):1T(p)=((4^(m+1-p))*T(p+1)-T(p))/(4^(m+1-p)-1); endm=m+1; endfprintf('数值积分结果为T(1)=%.8f\n',T(1)); fprintf('经过了%d 次迭代\n',m);用龙贝格积分法计算积分的近似值： 1.1请输入积分函数f=x^3请输入积分区间[a,b]=[6,100] 请输入误差限e=1000数值积分结果为T(1)=25000289.75465907 经过了15次迭代1.2请输入积分函数f=x^3请输入积分区间[a,b]=[6,100] 请输入误差限e=500数值积分结果为T(1)=24999982.87732918 经过了16次迭代2.1请输入积分函数f=sin(x)/x请输入积分区间[a,b]=[0.1e-010,1] 请输入误差限e=0.1e-004数值积分结果为T(1)=0.94607740 经过了12次迭代2.2请输入积分函数f=sin(x)/x请输入积分区间[a,b]=[0.1e-010,1] 请输入误差限e=0.1e-005数值积分结果为T(1)=0.94608236 经过了15次迭代w ww .k hd aw .c o m课后答案网3.1请输入积分函数f=sin(x^2) 请输入积分区间[a,b]=[0 1] 请输入误差限e=0.0001数值积分结果为T(1)=0.31031986 经过了11次迭代3.2请输入积分函数f=sin(x^2) 请输入积分区间[a,b]=[0,1] 请输入误差限e=0.1e-005数值积分结果为T(1)=0.31026911 经过了17次迭代实验六：常微分方程初值问题的数值解法程序1：R-K 和 修正Adams 预测-校正法的组合 syms f x y;f=input('请输入f(x,y)=');jx=input('请输入计算区间[a,b]='); yy(1)=input('请输入初值y0='); h=0.1;xx=jx(1):0.1:jx(2); for n=1:3x=xx(n); y=yy(n); k1=eval(f); x=xx(n)+h/2; y=yy(n)+k1*h/2; k2=eval(f);y=yy(n)+k2*h/2; k3=eval(f); x=xx(n+1); y=yy(n)+h*k3; k4=eval(f);yy(n+1)=yy(n)+(k1+2*k2+2*k3+k4)*h/6; endc(4)=yy(4); p(4)=yy(4); for i=1:4x=xx(i); y=yy(i); ff(i)=eval(f); endfor n=4:(length(xx)-1)p(n+1)=yy(n)+(55*ff(n)-59*ff(n-1)+37*ff(n-2)-9*ff(n-3))*h/24; m(n+1)=p(n+1)+251*(yy(n)-p(n))/270; x=xx(n+1); y=m(n+1);w ww .k hd aw .c o m课后答案网fff(n+1)=eval(f);c(n+1)=yy(n)+(9*fff(n+1)+19*ff(n)-5*ff(n-1)+ff(n-2))*h/24; yy(n+1)=c(n+1)-19*(c(n+1)-p(n+1))/270; y=yy(n+1); ff(n+1)=eval(f); endplot(xx,yy,'*') hold on程序2：Runge-Kutta 法 syms f x y;f=input('请输入f(x,y)=');jx=input('请输入计算区间[a,b]='); yy(1)=input('请输入初值y0='); h=0.1;xx=jx(1):0.1:jx(2); for n=1:(length(xx)-1) x=xx(n); y=yy(n); k1=eval(f); x=xx(n)+h/2; y=yy(n)+k1*h/2; k2=eval(f);y=yy(n)+k2*h/2; k3=eval(f); x=xx(n+1); y=yy(n)+h*k3; k4=eval(f);yy(n+1)=yy(n)+(k1+2*k2+2*k3+k4)*h/6; endplot(xx,yy,'o')程序3：第一题的求解析解的程序 syms sv x jx=[-1,0];sv=dsolve('Dy=(x^2-y^2)','y(-1)=0','x'); xxx=jx(1):0.01:jx(2); for n=1:length(xxx) x=xxx(n);yyy(n)=eval(sv); endplot(xxx,yyy)w ww .k hd aw .c o m课后答案网程序4：第二题的求解析解的程序 syms sv x jx=[0,1.5];sv=dsolve('Dy=(y-2*x/y)','y(0)=1','x'); xxx=jx(1):0.01:jx(2); for n=1:length(xxx) x=xxx(n);yyy(n)=eval(sv); endplot(xxx,yyy)图1：第一题的计算图形结果图2：第二题的计算图形结果：w ww .k hd aw .c o m课后答案网w ww .k hd aw .c o m课后答案网。
矩阵与数值分析上机实习题汇总
矩阵与数值分析上机实习题汇总矩阵与数值分析上机实习1.设, 其精确值为.(1)编制按从⼤到⼩的顺序, 计算的通⽤程序(2)编制按从⼩到⼤的顺序, 计算的通⽤程序(3)按两种顺序分别计算并指出有效位数（编制程序时⽤单精度）(4)通过本上机题，你明⽩了什么从⼩到⼤,代码:%1---SN = %N = input('please input a number(N>=2)')if(N < 2)disp('wrong number')elseS = 0;for j = 2:1:NS = S + 1/(j^2 -1);enddisp('S:')disp(S)end结果please input a number(N>=2)10^2N =100S:7.4005e-001>> clearplease input a number(N>=2)10^4N =10000S:7.4990e-001>> clearplease input a number(N>=2)10^6N =1000000S:7.5000e-001>>从⼤到⼩代码:%1---SN = %eps('single')N = input('please input a number(N>=2)') if(N < 2) disp('wrong number')elseS = 0;for j = N:-1:2S = S + 1/(j^2 -1);enddisp('S:')disp(S)end结果please input a number(N>=2)10^2N =100S:7.4005e-001>> clearans =1.1921e-007please input a number(N>=2)10^4N =10000S:7.4990e-001>> clearans =1.1921e-007please input a number(N>=2)10^6N =1000000S:7.5000e-001(4)计算的顺序影响结果。
数值分析实验答案
实验0截断误差与舍入误差#include ”stdio.h”#include ”math.h"const double 1112=0.693147190546;const double el=5e-6;void main(){mt sign;double s;long i；s=0.0;sign=l;i=l；wliile(fabs(hi2-s)>=el){s+=(1.0/i)*sign;sign=-sign:i++;}戸山1厅(”11=%1(1皿1);getch();}实验1拉格朗口插值法编写拉格朗口插值法通用子程序，并用以下函数表来上机求/(0.15)丿(0・31)。
#include <stdio.h>main(){static float Lx[10],Ly[10]; mt nj j;float xyp; pnntf(M enter n=M); scanfC%cT,&n); pnntf(M enter xfoT);fbr(i=0;i<n;i-H-) scanf(H%f\&Lx[i]);pnntf(M enter yi\n n); fbr(i=0;i<n;i-H-) scanf(H%f\&Ly[i]);pnntf(M enter x=”)；scanf(n%f\&x);n=6;Lx[O]=O;Lx[l]=0.1;Lx[2]=0.195;Lx[3]=0・3;Lx[4]=0.401;Lx[5]=0.5;Ly[0]=0.39894;Ly[l]=0.39695;Ly[2]=0.39142;Ly[3]=0・38138;Ly[4]=0.36812;Ly[5]=0.35206; x=0.15;*/fbr( i=O;i<n;i++){P=l；for(j=Oj<nj++){軟=J)p=p*(x-Lx[j])/(Lx[i]-Lx[j]);}y+=p*Ly[i];}prmtf(” y=%f\n",y);getch();}实验2最小二乘法测得铜导线在温度77(°C)时的电阻Ri0)如下表，求电阻与温度T的近似函数关系。
贵州大学数值分析往年试题(卷)(6套)
贵州大学2009级工程硕士研究生考试试卷数值分析注意事项：1.请考生按要求在下列横线内填写姓名、学号和年级专业。
2.请仔细阅读各种题目的回答要求，在规定的位置填写答案。
3.不要在试卷上乱写乱画，不要在装订线内填写无关的内容。
4．满分100分，考试时间120分钟。
专业学号姓名一、（12分）用牛顿迭代法求3220--=x x 在区间[1.5,2]内的一个近似根，要求31||10-+-<k k x x 。
二、（20分）已知()f x 的一组实验数据如下：（1）用三次插值公式求(1.28)f 的近似值；（2）用中心差商微分公式，求(1.5)'ƒ与求(2.0)'ƒ的近似值。
三、（20分）设方程组12312312335421537++=-+=--⎧⎪⎨⎩+=⎪x x x x x x x x x（1）用列主法求解方程组；（2）构造使G-S 方法收敛的迭代法，并取(0)(0,0,0)=T x，求方程组的二次迭代近似解根。
四、（16分）将积分区间2等分，分别用复化梯形公式与复化辛普森公式求21⎰x e dx的近似值。
五、（9分）设3211⎛⎫= ⎪--⎝⎭A，31⎛⎫= ⎪-⎝⎭x，求2||||x；谱半径()s A及条件数1()cond A。
六、（16分）取步长0.1=h ，用Euler 预报-校正公式求微分方程024|2='=--⎧⎨=⎩x y y x y 的解()y x 在x =0.1与x =0.2处的近似值(2)(0.1)y ，(2)(0.2)y 。
七、（7分）设A 为非奇异矩阵，0≠b ，%x 是=Ax b 的近似解，x 是=Ax b的解，证明1||||||||.()||||||||--≤%%b Ax x x cond A b x 。
贵州大学2010级工程硕士研究生考试试卷A数值分析注意事项：1.请考生按要求在下列横线内填写姓名、学号和年级专业。
2.请仔细阅读各种题目的回答要求，在规定的位置填写答案。
《数值分析》上机习题
《数值分析》上机习题1.用Newton法求方程X7 - 28X4 + 14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
#include<stdio.h>#include<math.h>int main(){ float x1,x,f1,f2;static int count=0;x1=0.1 ;//定义初始值do{x=x1;f1=x*(x*x*x*x*x*x-28*x*x*x)+14;f2=7*x*x*x*x*x*x-4*28*x*x*x;//对函数f1求导x1=x-f1/f2; count++; }while(fabs(x1-x)>=1e-5);printf("%8.7f\n",x1); printf("%d\n",count);return 0;}2.已知函数值如下表：试用三次样条插值求f(4.563)及f’(4.563)的近似值。
include <stdio.h>#include <math.h>#define N 11main(){double B[N+1][N+1],m,x,u[N+1],y[N+1],c[N+1],d[N+1];double e[N+1]={2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.675460 6,12.47667,13.1833476,13.8155106,14.0155106};int i;x=4.563;B[0][0]=-2;B[0][1]=-4;B[N][N-1]=4;B[N][N]=2;for(i=1;i<N;i++){B[i][i-1]=1;B[i][i]=4;B[i][i+1]=1;}u[0]=B[0][0];y[0]=e[0];for(i=1;i<N;i++){m=B[i][i-1]/u[i-1];u[i]=B[i][i]-m*B[i-1][i];y[i]=e[i]-m*y[i-1];}c[N]=y[N]/u[N];for(i=N-1;i>=0;i--)c[i]=(y[i]-B[i][i+1]*c[i+1])/u[i];for(i=0;i<12;i++){m=fabs(x-i);if(m>=2)d[i]=0;else if(m<=1)d[i]=0.5*fabs(pow(m,3))-m*m+2.0/3;elsed[i]=(-1.0/6.0)*fabs(pow(m,3))+m*m-2*fabs(m)+4/3.0;} m=0;for(i=0;i<12;i++) m=m+c[i]*d[i];printf("f(%4.3f)=%f\n",x,m);printf("f'(4.563)=%lf\n",(c[4]-c[2])/2); }3.用Romberg 算法求)00001.0(sin )75(32314.1=+⎰ε允许误差dx x x x x . #include "stdafx.h" #include<stdio.h> #include<math.h> float f(float x) {float f=0.0;f=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return (f); } main() {int i=1,j,k,n=12;float T[12],a=1.0,b=3.0,s=0.0; T[0]=0.5*(b-a)*(f(a)+f(b));for(j=1;j<n-1;j++){ for(k=1;k<=pow(2,j-1);k++) s+=f(a+(2*k-1)*(b-a)/pow(2,j)); T[j]=0.5*(T[j-1]+(b-a)*s/pow(2,j-1)); s=0.0; }T[11]=(4*T[1]-T[0])/(float)3;for(;fabs(T[11]-T[0])>0.00001;i++) {T[0]=T[11];for(j=1;j<n-1-i;j++) T[j]=(pow(4,i)*T[j+1]-T[j])/(pow(4,i)-1); T[11]=(pow(4,i+1)*T[1]-T[0])/(pow(4,i+1)-1); }printf("%f\n",T[11]); }4. 用定步长四阶Runge-Kutta 求解一、程序要求用定步长四阶法求解 y1’=1y2’=y3 y3’=1000-1000y2-100y3(y1(0)=y2(0)=y3(0)=0) h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3)h =0.0005，打印y i (0.025) ，⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧===--===0)0(0)0(0)0(10010001000//1/321323321y y y y y dt dy ydt dy dt dyy i(0.045) ，y i(0.085) ，y i(0.1) ，(i=1，2，3)#include "stdafx.h"#include <stdio.h>#include <math.h>double F(double x,double y[4],double f[4]){f[1]=0*x+0*y[1]+0*y[2]+0*y[3]+1;f[2]=0*x+0*y[1]+0*y[2]+1*y[3]+0;f[3]=0*x+0*y[1]-1000*y[2]-100*y[3]+1000;return(1);}void main(){double F(double x,double y[4],double f[4]);doubleh=0.0005,x=0,Y[4],k[5][4],s[4],f[4],det,m[4]={0.025,0.045,0.085,0.1};int i,j,t; for(t=0;t<=3;t++)/*龙格-库塔算法*/{for(j=0;j<=3;j++)Y[j]=0; //每求一组值后将初值条件还原为0 for(i=1;i<=int(m[t]/h);i++){for(j=1;j<=3;j++)s[j]=Y[j];det=F(x,s,f);for(j=1;j<=3;j++)k[1][j]=h*f[j]; /*四阶古典公式中的k值和求和的计算*/ for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[1][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[2][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[2][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[3][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+k[3][j];det=F(x+h,s,f);for(j=1;j<=3;j++)k[4][j]=h*f[j];for(j=1;j<=3;j++)Y[j]=Y[j]+(k[1][j]+2*k[2][j]+2*k[3][j]+k[4][j])/6;x+=h;} for(j=1;j<=3;j++)printf("y[%d](%f)=%f ",j,m[t],Y[j]); printf("\n"); } } 5.⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=40.00001 4.446782 2.213474- 0.238417 1.784317 0.037585- 1.010103- 3.123124 2.031743- 4.446782 30.719334 3.123789 1.103456- 2.121314 0.71828- 0.336993 1.112348 3.067813 2.213474- 3.123789 14.7138465 0.103458- 3.111223- 2.101023 1.056781- 0.784165- 1.7423820.238417 1.103456- 0.103458- 9.789365 0.431637 3.741856- 1.836742 1.563849 0.718719 1.7843172.1213143.111223- 0.431637 19.8979184.101011 2.031454 2.189736 0.113584-0.037585- 0.71828- 2.101023 3.741856- 4.101011 27.108437 3.123848 1.012345- 1.112336 1.010103- 0.336993 1.056781- 1.836742 2.031454 3.123848 15.567914 3.125432- 1.061074- 3.123124 1.112348 0.784165- 1.563849 2.189736 1.012345- 3.125432- 19.141823 2.115237 2.031743- 3.067813 1.742382 0.718719 0.113584- 1.112336 1.061074- 2.115237 12.38412A Tb )5.6784392- 4.719345 1.1101230 86.612343- 1.784317 0.84671695 25.173417- 33.992318 2.1874369(= 用列主元消去法求解Ax=b 。
贵州大学数值分析往年试题(6套)
贵州大学2020级工程硕士研究生考试试卷数值分析注意事项：1.请考生按要求在以下横线内填写姓名、学号和年级专业。
2.请认真阅读各类题目的回答要求，在规定的位置填写答案。
3.不要在试卷上乱写乱画，不要在装订线内填写无关的内容。
4．总分值100分，考试时刻120分钟。
专业 学号 姓名一、（12分）用牛顿迭代法求3220--=x x 在区间[1.5,2]内的一个近似根，要求31||10-+-<k k x x 。
（1）用三次插值公式求(1.28)f 的近似值；（2）用中心差商微分公式，求(1.5)'ƒ与求(2.0)'ƒ的近似值。
三、（20分）设方程组12312312335421537++=-+=--⎧⎪⎨⎩+=⎪x x x x x x x x x（1）用列主法求解方程组；（2）构造使G-S 方式收敛的迭代法，并取(0)(0,0,0)=T x，求方程组的二次迭代近似解根。
四、（16分）将积分区间2等分，别离用复化梯形公式与复化辛普森公式求21⎰x e dx的近似值。
五、（9分）设3211⎛⎫= ⎪--⎝⎭A，31⎛⎫= ⎪-⎝⎭x，求2||||x；谱半径()s A及条件数1()cond A。
六、（16分）取步长0.1=h ，用Euler 预报-校正公式求微分方程024|2='=--⎧⎨=⎩x y y x y 的解()y x 在x =0.1与x =0.2处的近似值(2)(0.1)y ，(2)(0.2)y 。
七、（7分）设A 为非奇异矩阵，0≠b ，x 是=Ax b 的近似解，x 是=Ax b的解，证明1||||||||.()||||||||--≤b Ax x x cond A b x 。
贵州大学2020级工程硕士研究生考试试卷A数值分析注意事项：1.请考生按要求在以下横线内填写姓名、学号和年级专业。
2.请认真阅读各类题目的回答要求，在规定的位置填写答案。
3.不要在试卷上乱写乱画，不要在装订线内填写无关的内容，4．总分值100分，考试时刻120分钟。
数值分析上机实习题
数值分析上机实习题第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值画出插值函数及图形。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析上机实验报告课程名称：数值分析上机实验学院：机械工程学院专业：机械制造姓名：******* 学号：********** 年级：12级任课教师：***老师2012年12月30日一．已知A 与b12.38412 2.115237 -1.061074 1.112336 -0.1135840.718719 1.742382 3.067813 -2.031743 2.11523719.141823 -3.125432 -1.012345 2.189736 1.563849-0.784165 1.112348 3.123124 -1.061074 -3.125A =43215.567914 3.123848 2.031454 1.836742-1.056781 0.336993 -1.010103 1.112336 -1.012345 3.12384827.108437 4.101011-3.741856 2.101023 -0.71828 -0.037585 -0.1135842.189736 2.031454 4.10101119.8979180.431637-3.111223 2.121314 1.784317 0.718719 1.563849 1.836742 -3.741856 0.4316379.789365-0.103458 -1.103456 0.238417 1.742382 -0.784165 -1.056781 2.101023-3.111223-0.10345814.7138465 3.123789 -2.213474 3.067813 1.112348 0.336993-0.71828 2.121314-1.103456 3.12378930.719334 4.446782 -2.031743 3.123124-1.010103-0.0375851.7843170.238417-2.213474 4.44678240.00001[ 2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.11012304.719345 -5.6784392]T B ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦=（2）用超松弛法求解Bx=b （取松弛因子ω=1.4，x （0）=0，迭代9次）。
（3）用列主元素消去法求解Bx=b 。
解：（3）、用列主元素消去法求解Bx=b（一）、理论依据：其基本思想是选取绝对值尽量大的元素作为主元素，进行行与列的交换，再进行回代，求出方程的解。
将方阵A 和向量b 写成C=（A b ）。
将C 的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素1j c ，将第j 行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n 个元素都消成0。
将变换后的矩阵(1)C 的第二列中第二行的元素与其下面的此列的元素逐一进行比较，找到最大的元素(1)2k c ，将第k 行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n 个元素都消成0。
以此方法将矩阵的左下部分全都消成0。
（二）、计算程序： #include "math.h" #include "stdio.h" void main() { double u[9],x1[9],y[9],q[9],b1[9][10],x[9],a[9][9]={{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103}, {1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585}, {-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317}, {0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417}, {1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474}, {3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}};int sign(double x);double k,t,s,w,e,c,z;int i,j,n,r;doubleb[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};for(r=0;r<=6;r++) /*Household 变换*/{e=0.0;for(i=r+1;i<=8;i++)e=e+a[i][r]*a[i][r];s=sqrt(e);t=s*s+fabs(a[r+1][r])*s;for(i=0;i<=r;i++) u[i]=0;c=a[r+1][r]; /*求u[i]的值*/u[r+1]=a[r+1][r]+s*sign(c);for(i=r+2;i<=8;i++)u[i]=a[i][r];for(i=0;i<=8;i++){y[i]=0;for(j=0;j<=8;j++)y[i]+=a[i][j]*u[j]/t;} /*求出y向量*/ k=0;for(i=0;i<9;++i)k+=0.5*(u[i]*y[i])/t;for(i=0;i<=8;i++)q[i]=y[i]-k*u[i];for(i=0;i<=8;i++)for(j=0;j<=8;j++)a[i][j]=a[i][j]-(q[j]*u[i]+u[j]*q[i]);} /*求结果*/printf("Household变换:\n");for(i=0;i<9;++i)for(j=0;j<9;++j) /*打印转化后的矩阵*/ {if (j%9==0)printf("\n");printf("%-9.5f",a[i][j]);}printf("\n");printf("超松弛变量法得解:\n");w=1.4; /*超松弛法*/for(i=0;i<9;i++)x1[i]=0;for(i=0;i<9;i++)for(j=0;j<9;j++){if(i==j)b1[i][j]=0;else b1[i][j]=-a[i][j]/a[i][i];} /*求出矩阵b1[9][9]和b1[i][9]的值*/ for(i=0;i<9;i++)b1[i][9]=b[i]/a[i][i];for(n=0;n<9;n++)for(i=0;i<9;i++){z=0;for(j=0;j<9;j++)z=z+b1[i][j]*x1[j]; /*执行本算法*/z=z+b1[i][9];x1[i]=x1[i]*(1-w)+w*z;}for(i=0;i<9;i++){if (i==5)printf("\n");printf("x%d=%-10.6f",i,x1[i]);}printf("\n");printf("列主元消去法得解:\n");u[0]=a[0][0]; /*以下是消去法*/y[0]=b[0];for(i=1;i<9;i++){q[i]=a[i][i-1]/u[i-1];u[i]=a[i][i]-q[i]*a[i-1][i];y[i]=b[i]-q[i]*y[i-1];}x[8]=y[8]/u[8]; /*执行本算法*/for(i=7;i>=0;i--)x[i]=(y[i]-a[i][i+1]*x[i+1])/u[i]; /*求出x的值*/ for(i=0;i<9;i++){if (i==5)printf("\n");printf(" x%d=%-10.6f",i,x[i]);}printf("\n");}int sign(double x){int z;z=(x>=(1e-6)?1:-1);return(z);}（三）、计算结果打印：（四）、问题讨论：a ．由于选主元，使|lij|最小，这样去乘方程的每一系数时，系数的舍入误差不至扩大，并防止溢出与停机。
b ． 由于选主元，回代时作除数主元aii(i)的绝对值也是最大,这样扩大的误差也是最小.c ． 若detA ≠0,则选主元后的主元aii(i) ≠0,这是因为若aii(i)=0则必有aji(i)=0 (j>i)这样按行列式的Laplace 展开式就有detA=0而矛盾,因此用主元消去法中进行下去,不止中断停机.同时也由于选主元, aii(i)接近零的概率减少.运行结果基本与准确值无异，因为这种算法的无误差的算法。
三．试用三次样条插值求()4.563f 及(4.563)f '的近似值。
（一）方法由s(x i )=y i ,I=1,2,…N s’(x 0)=y 0’,s ’(x N )=y N ’可得S(x i )= ∑+-=11N j c j Ω3(x i -x j /h)=y iS ’(x 0)=1/h ∑+-=11N j c j Ω3’(x 0-x j /h)=y ’0S ’(x N )=1/h ∑+-=11N j c j Ω3’(x N -x j /h)=y ’N其中j c 必须的方程 ：⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛---2.081551.1318338.1347667.12675460.11750557.106566268.93177664.85916738.6158839.4021011411411411411411411411411411411011098765432101c c c c c c c c c c c cq[0]=0 , u[0]=0 ,1,,2,1]),[][][(][][-⋅⋅⋅=⋅+-=n i i q i a i b i c i qn i i q i a i b i u i a i d i u ,,2,1]),1[][][(])[][][(}[⋅⋅⋅=-⋅+⋅-= x[9]=u[9]1,,2,1],[]1[][][⋅⋅⋅--=++⋅=n n i i u i x i q i x []333333)2()1(46)1(4)2(61)(+++++-+--++-+=Ωx x x x x x s ' (x)= ∑∑-=-=--Ω--+Ω10121012)21()21(j j j j j x c j x c（二）原程序及运行结果：#include<stdio.h> #include<math.h> #define N 9 void main() {double x[N+1]={1,2,3,4,5,6,7,8,9,10},y[N+1]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1 972246,2.3025851},h[N+1],d[N+1],a[N+1],c[N+1],b[N+1]={2,2,2,2,2,2,2,2,2,2},s[N+1],t[N+1],l[N+1],M[N+1], f,f1;int i;for(i=1;i<=N;i++) /*使用对任意分化的三弯矩插值法*/h[i-1]=x[i]-x[i-1];d[0]=6/h[0]*((y[1]-y[0])/h[0]-1);d[N]=6/h[N-1]*(0.1-(y[N]-y[N-1])/h[N-1]);for(i=1;i<=N-1;i++){d[i]=6/(h[i-1]+h[i])*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]);a[i]=h[i-1]/(h[i-1]+h[i]);c[i]=1-a[i];}c[0]=1;a[N]=1;s[0]=b[0];t[0]=d[0]; /*用追赶法求解三对角方程组*/for(i=0;i<=N-1;i++){l[i+1]=a[i+1]/s[i];s[i+1]=b[i+1]-l[i+1]*c[i];t[i+1]=d[i+1]-l[i+1]*t[i];}M[N]=t[N]/s[N];for(i=N-1;i>=0;i--)M[i]=(t[i]-c[i]*M[i+1])/s[i];f=M[3]*pow((x[4]-4.563),3)/6/h[3] /*求算4.563这点的函数值*/+M[4]*pow((4.563-x[3]),3)/6/h[3]+(y[3]-M[3]*h[3]*h[3]/6)*(x[4]-4.563)/h[3]+(y[4]-M[4]*h[3]*h[3]/6)*(4.563-x[3])/h[3];f1=-3*M[3]*pow((x[4]-4.563),2)/6/h[3] /*求算4.563这点的一阶导数值*/+3*M[4]*pow((4.563-x[3]),2)/6/h[3] -(y[3]-M[3]*h[3]*h[3]/6)/h[3] +(y[4]-M[4]*h[3]*h[3]/6)/h[3];printf("f(4.563)=%lf f'(4.563)=%lf\n",f,f1);}（三） 输出结果：（四） 问题讨论：其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中不取1,x,,x 2,x 3,(x-x j )+3为基函数,而取B 样条函数Ω3(x-x j /h)为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X -1,X N+1,则任意三次样条函数可用Ω3(x-x j /h)线性组合来表示 S(x)= ∑+-=11N j c j Ω3(x-x j /h) 这样对不同插值问题,若能确定c j 由解的唯一性就能求得解S(x).同时样条插值效果比Lagrange 插值好,近似误差较小.没有Runge 现象。