《数值分析》计算实习题
姓名：
学号：
班级：
第二章
1、程序代码
Clear;clc;
x1=[0.2
0.4
0.6
0.8
1.0];
y1=[0.98
0.92
0.81
0.64
0.38];
n=length(y1);
c=y1(:);
for
j=2:n
%求差商
for
i=n:-1:j
c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));
end
end
syms
x
df
d;
df(1)=1;d(1)=y1(1);
for
i=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=pp.coefs;
q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1];
q1=vpa(collect(q1),5)
q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1];
q2=vpa(collect(q2),5)
q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1];
q3=vpa(collect(q3),5)
q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1];
q4=vpa(collect(q4),5)%求解并化简多项式
2、运行结果
P4
=
0.98*x
0.3*(x
0.2)*(x
0.4)
0.625*(x
0.2)*(x
0.4)*(x
0.6)
0.20833*(x
0.2)*(x
0.4)*(x
0.8)*(x
0.6)
+
0.784
q1
=
1.3393*x^3
+
0.80357*x^2
0.40714*x
+
1.04
q2
=
1.3393*x^3
+
1.6071*x^2
0.88929*x
+
1.1643
q3
=
1.3393*x^3
+
2.4107*x^2
1.6929*x
+
1.4171
q4
=
1.3393*x^3
+
3.2143*x^2
2.8179*x
+
1.86293、问题结果
4次牛顿差值多项式=
0.98*x
0.3*(x
0.2)*(x
0.4)
0.625*(x
0.2)*(x
0.4)*(x
0.6)
0.20833*(x
0.2)*(x
0.4)*(x
0.8)*(x
0.6)
+
0.784
三次样条差值多项式
第三章
1、程序代码
Clear;clc;
x=[0
0.1
0.2
0.3
0.5
0.8
1];
y=[1
0.41
0.5
0.61
0.91
2.02
2.46];
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
=
-6.6221
12.8147
-4.6591
0.9266
p2
=
2.8853
-12.3348
16.2747
-5.2987
0.9427
p3
=
3.1316
-1.2400
0.73563、问题结果
三次多项式拟合P1=
四次多项式拟合P2=
二次多项式拟合P3=
第四章
1、程序代码
1)建立函数文件f.m:
function
y=fun(x);
y=sqrt(x)*log(x);
2)编写程序：
a.利用复化梯形公式及复化辛普森公式求解：
Clear;clc;
h=0.001;%h为步长,可分别令h=1,0.1,0.01,0.001
n=1/h;t=0;s1=0;s2=0;
for
i=1:n-1
t=t+f(i*h);
end
T=h/2*(0+2*t+f(1));T=vpa(T,7)
%梯形公式
for
i=0:n-1
s1=s1+f(h/2+i*h);
end
for
i=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:0.001:1;
%h为步长,可分别令h=1,0.1,0.01,0.001
y=sqrt(x).*log(x+eps);
T=trapz(x,y);
T=vpa(T,7)
（只是h=1的运行结果不一样，T=1.110223*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;
for
i=1:m
h=h/2;n=1/h;t=0;
for
j=1:n-1
t=t+f(j*h);
end
T(i)=h/2*(0+2*t+f(1));%梯形公式
end
for
i=1:m-1
for
j=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
=
-0.4443875
S
=
-0.4444345
ans
=
-0.44444143、问题结果
a.利用复合梯形公式及复合辛普森公式求解：
步长h
0.1
0.01
0.001
梯形求积T=
0
[1.110223*10^(-16)]
-0.4170628
-0.4431179
-0.4443875
辛普森求积S=
-0.3267527
-0.4386308
-0.4441945
-0.4444345
b.利用龙贝格公式求解：
通过15次二分，得到结果：-0.4444414
第五章
1、程序代码
（1）LU分解解线性方程组：
Clear;clc;
A=[10
0
2.099999
0
2];
b=[8;5.900001;5;1];
[m,n]=size(A);
L=eye(n);
U=zeros(n);
flag=＇ok＇;
for
i=1:n
U(1,i)=A(1,i);
end
for
r=2:n
L(r,1)=A(r,1)/U(1,1);
end
for
i=2:n
for
j=i:n
z=0;
for
r=1:i-1
z=z+L(i,r)*U(r,j);
end
U(i,j)=A(i,j)-z;
end
if
abs(U(i,i)) flag=＇failure＇ return; end for k=i+1:n m=0; for q=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）列主元消去法： function x = gauss(A,b); A=[10 0 1;-3 2.099999 2;5 -1;2 0 2]; b=[8;5.900001;5;1]; [n,n] = size(A); x = zeros(n,1); Aug = [A,b]; %增广矩阵 for k = 1:n-1 [piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r r = r + k 1; % 列主元所在大矩阵的行 if r>k temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)==0,error(‘对角元出现0’),end % 把增广矩阵消元成为上三角 for p = 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); for k = n-1:-1:1 x(k)=b(k); for p=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 = 1.0e+006 * 0.0000 0 0 0 -0.0000 0.0000 0 0 0.0000 -2.5000 0.0000 0 0.0000 -2.4000 0.0000 0.0000 U = 1.0e+007 * 0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 0 0 0.0000 x = -0.0000 -1.0000 1.0000 1.0000 detA = -762.0001 2）列主元消去法 detA = 762.0001 ans = 0.0000 -1.0000 1.0000 1.00003、问题结果 1） LU分解解线性方程组 L= U= x=（0.0000，-1.0000，1.0000，1.0000）T detA=-762.001 2)列主元消去法 x=（0.0000，-1.0000，1.0000，1.0000）T detA=762.001 第六章 1、程序代码 （1）Jacobi迭代 Clear;clc; n = 6; %也可取n=8,10 H = hilb(n); b = H * ones(n,1); e = 0.00001; for i = 1:n if H(i,i)==0 ＇对角元为零，不能求解＇ return end end x = zeros(n,1); k = 0; kend = 10000; r = 1; while k<=kend & r>e x0 = x; for i = 1:n s = 0; for j = 1:i s = s + H(i,j) * x0(j); end for j = 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 if k>kend ＇迭代不收敛，失败＇ else ＇求解成功＇ x k end （2）SOR迭代 1）程序代码 function s = SOR(n,w); H = hilb(n); b = H*ones(n,1); e = 0.00001; for i = 1:n if H(i,i)==0 ‘对角线为零，不能求解’ return end end x = zeros(n,1); k = 0; kend = 10000; r = 1; while k<=kend & r>e x0 = x; for i = 1:n s = 0; for j = 1:i s = s + H(i,j) * x(j); end for j = 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 if k>kend ＇迭代不收敛，失败＇ else ＇求解成功＇ x end 2） 从命令窗口中分别输入： SOR（6,1） SOR（8,1） SOR（10,1） SOR（6,1.5） SOR（8,1.5） SOR（10，1.5） 2、运行结果 Jacobi迭代： ans = 迭代不收敛，失败 SOR迭代： 第七章 1、程序代码 (1)不动点迭代法 1）建立函数文件：g.m function f=g(x) f(1)=20/(x^2+2*x+10); 2)建立函数文件：bdd.m function [y,n] = bdd(x,eps) if nargin==1 eps=1.0e-8; elseif nargin<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; format long; 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=0.5; disp(x1); Eps=1E-8; k=0; while x0=x1; k=k+1; x1=feval(Fx,x1); %将x1代入牛顿迭代公式替代x1 disp(x1); %在屏幕上显示x1 if k==m break; end end k,x12、运行结果 （1）不动点迭代法 >> bdd(0) n = ans = 1.3688 （2） 牛顿迭代 x= 0 1.*** 1.37***21 1.*** 1.*** 1.*** 1.*** 1.*** k = x1 = 1.*** 3、问题结果 (1)不动点迭代法 x=1.3688 n=25 收敛太慢 （2）牛顿迭代 初值取0 迭代次数k=8时，k x(k) 1.4666666 1.3715120 1.3688102 1.3688081 1.3688081 1.3688081 1.3688081
点击咨询