数值分析计算实习题.docx
- 文档编号:11661427
- 上传时间:2023-03-29
- 格式:DOCX
- 页数:22
- 大小:116.96KB
数值分析计算实习题.docx
《数值分析计算实习题.docx》由会员分享,可在线阅读,更多相关《数值分析计算实习题.docx(22页珍藏版)》请在冰豆网上搜索。
数值分析计算实习题
-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN
数值分析计算实习题
《数值分析》计算实习题
姓名:
学号:
班级:
第二章
1、程序代码
Clear;clc;
x1=[];
y1=[];
n=length(y1);
c=y1(:
);
forj=2:
n%求差商
fori=n:
-1:
j
c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));
end
end
symsxdfd;
df
(1)=1;d
(1)=y1
(1);
fori=2:
n%求牛顿差值多项式
df(i)=df(i-1)*(x-x1(i-1));
d(i)=c(i-1)*df(i);
end
P4=vpa(sum(d),5)%P4即为4次牛顿插值多项式,并保留小数点后5位数
pp=csape(x1,y1,'variational');%调用三次样条函数
q=;
q1=q(1,:
)*[^3;^2;;1];
q1=vpa(collect(q1),5)
q2=q(1,:
)*[^3;^2;;1];
q2=vpa(collect(q2),5)
q3=q(1,:
)*[^3;^2;;1];
q3=vpa(collect(q3),5)
q4=q(1,:
)*[^3;^2;;1];
q4=vpa(collect(q4),5)%求解并化简多项式
2、运行结果
P4=
*x-*(x-*(x--*(x-*(x-*(x--*(x-*(x-*(x-*(x-+
q1=
-*x^3+*x^2-*x+
q2=
-*x^3+*x^2-*x+
q3=
-*x^3+*x^2-*x+
q4=
-*x^3+*x^2-*x+
3、问题结果
4次牛顿差值多项式
=*x-*(x-*(x--*(x-*(x-*(x--*(x-*(x-*(x-*(x-+
三次样条差值多项式
第三章
1、程序代码
Clear;clc;
x=[01];
y=[1];
p1=polyfit(x,y,3)%三次多项式拟合
p2=polyfit(x,y,4)%四次多项式拟合
y1=polyval(p1,x);
y2=polyval(p2,x);%多项式求值
plot(x,y,'c--',x,y1,'r:
',x,y2,'y-.')
p3=polyfit(x,y,2)%观察图像,类似抛物线,故用二次多项式拟合。
y3=polyval(p3,x);
plot(x,y,'c--',x,y1,'r:
',x,y2,'y-.',x,y3,'k--')%画出四种拟合曲线
2、运行结果
p1=
p2=
p3=
3、问题结果
三次多项式拟合P1=
四次多项式拟合P2=
二次多项式拟合P3=
第四章
1、程序代码
1)建立函数文件:
functiony=fun(x);
y=sqrt(x)*log(x);
2)编写程序:
a.利用复化梯形公式及复化辛普森公式求解:
Clear;clc;
h=;%h为步长,可分别令h=1,,,
n=1/h;t=0;s1=0;s2=0;
fori=1:
n-1
t=t+f(i*h);
end
T=h/2*(0+2*t+f
(1));T=vpa(T,7)%梯形公式
fori=0:
n-1
s1=s1+f(h/2+i*h);
end
fori=1:
n-1
s2=s2+f(i*h);
end
S=h/6*(0+4*s1+2*s2+f
(1));S=vpa(S,7)%辛普森公式
a’复化梯形公式和复化辛普生公式程序代码也可以是:
Clear;clc;
x=0:
:
1;%h为步长,可分别令h=1,,,
y=sqrt(x).*log(x+eps);
T=trapz(x,y);
T=vpa(T,7)
(只是h=1的运行结果不一样,T=*10^(-16),而其余情况结果都相同)
Clear;clc;
f=inline('sqrt(x).*log(x)',x);
S=quadl(f,0,1);
S=vpa(S,7)
b.利用龙贝格公式求解:
Clear;clc;
m=14;%m+1即为二分次数
h=2;
fori=1:
m
h=h/2;n=1/h;t=0;
forj=1:
n-1
t=t+f(j*h);
end
T(i)=h/2*(0+2*t+f
(1));%梯形公式
end
fori=1:
m-1
forj=m:
i+1
T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);
%通过不断的迭代求得T(j),即T表的对角线上的元素。
end
end
vpa(T(m),7)
2、运行结果
T=
S=
ans=
3、问题结果
a.利用复合梯形公式及复合辛普森公式求解:
步长h
1
梯形求积T=
0
[*10^(-16)]
辛普森求积S=
b.利用龙贝格公式求解:
通过15次二分,得到结果:
第五章
1、程序代码
(1)LU分解解线性方程组:
Clear;clc;
A=[10-701
-362
5-15-1
2102];
b=[8;;5;1];
[m,n]=size(A);
L=eye(n);
U=zeros(n);
flag='ok';
fori=1:
n
U(1,i)=A(1,i);
end
forr=2:
n
L(r,1)=A(r,1)/U(1,1);
end
fori=2:
n
forj=i:
n
z=0;
forr=1:
i-1
z=z+L(i,r)*U(r,j);
end
U(i,j)=A(i,j)-z;
end
ifabs(U(i,i)) flag='failure' return; end fork=i+1: n m=0; forq=1: i-1 m=m+L(k,q)*U(q,i); end L(k,i)=(A(k,i)-m)/U(i,i); end end L U y=L\b;x=U\y detA=det(L*U) (2)列主元消去法: functionx=gauss(A,b); A=[10-701;-362;5-15-1;2102]; b=[8;;5;1]; [n,n]=size(A); x=zeros(n,1); Aug=[A,b];%增广矩阵 fork=1: n-1 [piv,r]=max(abs(Aug(k: n,k)));%找列主元所在子矩阵的行r r=r+k-1;%列主元所在大矩阵的行 ifr>k temp=Aug(k,: ); Aug(k,: )=Aug(r,: ); Aug(r,: )=temp; end ifAug(k,k)==0,error(‘对角元出现0’),end %把增广矩阵消元成为上三角 forp=k+1: n Aug(p,: )=Aug(p,: )-Aug(k,: )*Aug(p,k)/Aug(k,k); end end %解上三角方程组 A=Aug(: 1: n);b=Aug(: n+1); x(n)=b(n)/A(n,n); fork=n-1: -1: 1 x(k)=b(k); forp=n: -1: k+1 x(k)=x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end detA=det(A) 2、运行结果 1)LU分解解线性方程组 L= +006* 000 00 0 U= +007* 0 0 00 000 x= detA= 2)列主元消去法 detA= ans= 3、问题结果 1)LU分解解线性方程组 L= U= x= detA= 2)列主元消去法 x= detA= 第六章 1、程序代码 (1)Jacobi迭代 Clear;clc; n=6;%也可取n=8,10 H=hilb(n); b=H*ones(n,1); e=; fori=1: n ifH(i,i)==0'对角元为零,不能求解' return end end x=zeros(n,1); k=0; kend=10000; r=1; whilek<=kend&r>e x0=x; fori=1: n s=0; forj=1: i-1 s=s+H(i,j)*x0(j); end forj=i+1: n s=s+H(i,j)*x0(j); end x(i)=b(i)/H(i,i)-s/H(i,i); end r=norm(x-x0,inf); k=k+1; end ifk>kend'迭代不收敛,失败' else'求解成功' x k end (2)SOR迭代 1)程序代码 functions=SOR(n,w); H=hilb(n); b=H*ones(n,1); e=; fori=1: n ifH(i,i)==0‘对角线为零,不能求解’ return end end x=zeros(n,1); k=0; kend=10000; r=1; whilek<=kend&r>e x0=x; fori=1: n s=0; forj=1: i-1 s=s+H(i,j)*x(j); end forj=i+1: n s=s+H(i,j)*x0(j); end x(i)=(1-w)*x0(i)+w/H(i,i)*(b(i)-s); end r=norm(x-x0,inf); k=k+1; end ifk>kend'迭代不收敛,失败' else'求解成功' x end 2) 从命令窗口中分别输入: SOR(6,1) SOR(8,1) SOR(10,1) SOR(6,) SOR(8,) SOR(10,) 2、运行结果 Jacobi迭代: ans= 迭代不收敛,失败 SOR迭代: 第七章 1、程序代码 (1)不动点迭代法 1)建立函数文件: functionf=g(x) f (1)=20/(x^2+2*x+10); 2)建立函数文件: function[y,n]=bdd(x,eps) ifnargin==1 eps=; elseifnargin<1 error return end x1=g(x); n=1; while(norm(x1-x)>=eps)&&(n<=10000) x=x1; x1=g(x); n=n+1; end y=x; n 3)从命令窗口输入: bdd(0) (2)牛顿迭代 clear;clc; formatlong; m=8;%m为迭代次数,可分别令m=2,4,6,8,10 x=sym('x'); f=sym('x^3+2*x^2+10*x-20'); df=diff(f,x); FX=x-f/df;%牛顿迭代公式 Fx=inline(FX); disp('x='); x1=; disp(x1); Eps=1E-8; k=0; while1 x0=x1; k=k+1; x1=feval(Fx,x1);%将x1代入牛顿迭代公式替代x1 disp(x1);%在屏幕上显示x1 ifk==m break; end end k,x1 2、运行结果 (1)不动点迭代法 >>bdd(0) n= 25 ans= (2)牛顿迭代 x= 0 2 k= 8 x1= 3、问题结果 (1)不动点迭代法 x=n=25收敛太慢 (2)牛顿迭代 初值取0 迭代次数k=8时, k x(k) 1 2 2 3 4 5 6 7 8
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 计算 实习