贵州大学数值分析上机实验答案
(完整版)哈工大-数值分析上机实验报告
（完整版）哈工大-数值分析上机实验报告实验报告题目： 非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法 分法和 Newton 法及改进的 Newton 法。
刖言：理:对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法：二分法和 Newton 法。
对于二分法，其数学实质就是说对于给定的待求解的方程 f(x)，其在［a,b ］上连续，f(a)f(b)5e-6)R=b-a;%求出误差 k=k+1; end x=c%给出解Newton 法及改进的Newton 法源程序: clear %%%%输入函数 f=input （'请输入需要求解函数>>','s'） %%% 求解f （x ）的导数 df=diff(f);掌握二分法与 Newton 法的基本原理和应用。
数学原c=(a+b)/2; if f12(a)*f12(c)>0; a=c; elseb=c; end%%改进常数或重根数miu=2; %%初始值x0x0=in pu t('i nput initial value xO>>'); k=0;%max=100;%最大迭代次数x0时否就是解while (abs(R)>1e-8)x1=xO-miu*eval(subs(f,'xO','x'))/eval(subs(df,'xO',if (eval(subs(f,'x0','x'))max;%给定值，认为迭代不收敛，重新输入初值x0,y/n>>','s');if strc mp (ss,'y')k;%给出迭代次数x=x0 ；%给出解结果分析和讨论:1.用二分法计算方程在［1，2］内的根。
数值分析上机作业(总)
数值分析上机实验一、解线性方程组直接法（教材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公式需要在全体积分限内采用较小间距才满足条件。
数值分析上机第七章+第九章
数值分析第三次上机实验报告学院班级：学生学号：学生姓名：同作者:实验日期：1.实验题目： P232 3.(1) 一、实验目的：设f(x)=1/x,(1)求f(x)在[1,2] 上的零次和一次最佳一致逼近多项式。
(2)求f(x)在[1,2] 上的零次和一次最佳平凡逼近多项式。
二、实验环境：1.matlab2014b/macOS Seirra2.G 楼机房三、实验内容及实验原理：1.零次最佳逼近多项式 原理1: ()()02M m P x +=所以f(x)=1/x 在[1,2]上的零次最佳一致逼近多项式()01132P 24x +== 原理2：()()()0000,,f P x ϕϕϕ=()101P x a a x =+f(x)=1/x 在[1,2]上的零次最佳平方逼近多项式()()()210020011,ln 2,dx f x P x dxϕϕϕ===⎰⎰ 2. 一次最佳逼近多项式 (1)一次最佳一致逼近多项式： 解：21'()f x x =- ，32''()0f x x =>∴ 1，2为交错点，设101P ()x a a x =+111()()12212f b f a a b a --===---且由112111'(),2f x x x =-=-=1111(1()()()322224f a f x a a xa-+++ =-==故得131P()42x x+=-(2)一次最佳平方逼近多项式解：设10101P(),1,x c c x xϕϕ=+==001000011111(,)(,)(,)=(,)(,)(,)c fc fϕϕϕϕϕϕϕϕϕϕ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦由0111,,x fxϕϕ===得：2001(,)1dxϕϕ==⎰221117(,)3x dxϕϕ==⎰2100113(,)(,)2xdxϕϕϕϕ===⎰2011(,)ln2f dxxϕ==⎰211(,)1f dxϕ==⎰得到法方程组：01013ln2237123c cc c⎧+=⎪⎪⎨⎪+=⎪⎩解之得：**01c8.9782,c0.4766==1P()8.97820.4766x x∴=+四、实验结果及其分析:经拟合结果无误。
数值分析实验报告 课后答案_【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课后答案网。
数值分析实验答案
实验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的近似函数关系。
数值分析上机实验指导书
“数值计算方法”上机实验指导书实验一 误差分析实验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的全部根；而函数 p o l y (vb = 的输出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）中的ε。
实验要求：（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析上机实验报告课程名称：数值分析上机实验学院：机械工程学院专业：机械制造姓名：******* 学号：********** 年级：12级任课教师：***老师2012年12月30日一．已知A 与b12.38412 2.115237 -1.061074 1.112336 -0.1135840.7187191.7423823.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.5679143.123848 2.0314541.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.1101230 4.719345 -5.6784392]TB ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦=（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,,x2,x3,(x-xj)+3为基函数,而取B样条函数Ω3(x-xj/h)为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用Ω3(x-xj/h)线性组合来表示 S(x)= ∑+-=11NjcjΩ3(x-xj/h) 这样对不同插值问题,若能确定cj由解的唯一性就能求得解S(x).同时样条插值效果比Lagrange插值好,近似误差较小.没有Runge现象。