矩阵与数值分析例题matlab仿真

解答：
应用原有迭代公式： 程序： clc clear x(1)=1; i=1; format while (f(x(i))~=0)
i=i+1; x(i)=sqrt(10/(x(i-1)+4)); end fprintf('方程在[1,1.5]内的解为：x=%f 迭代次数为：k=%d',x(i),i)
《矩阵与数值分析》 实验报告
第一题
题目：
一、对于数列1, 1 , 1 , 1 , 1 ,，有如下两种生成方式 3 9 27 81
1、首项为 a0
1，递推公式为 an

1 3
an1, n
1, 2, ；
2、前两项为 a0
1, a1

1 3
，递推公式为
an

10 3
an1
an2 , n
x0=x; x=B*x0+f; k=k+1; end fprintf('方程组的近似解为：') x fprintf('迭代次数为：%d 次',k) 运行结果：
Gauss-Seidel 迭代法： 程序： clc clear A=[6 2 1 -2;2 5 0 -2;-2 0 8 5;1 3 2 7]; b=[4;7;-1;0]; D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=inv(D-L)*U; f=inv(D-L)*b; x0=[0;0;0;0]; x=B*x0+f; k=0; while(max(abs(x(1:3)-x0(1:3)))>=1e-6)
a=abs(A(i:n,i)); d=max(a); l=find(a(1:(n-i+1))==d); l=l+i-1; c=A(l,:); A(l,:)=A(i,:); A(i,:)=c; for j=i+1:n
l1=A(j,i)/A(i,i); A(j,:)=A(j,:)-l1*A(i,:); end end x=zeros(n,1); x(n)=A(n,n+1)/A(n,n); for i=1:n-1 z=A(n-i,n+1); for j=n-i+1:n; z=z-A(n-i,j)*x(j); end x(n-i)=z/A(n-i,n-i); end fprintf(‘方程组的解为：\n’); x 运行结果：
end f=fval; subplot(2,3,k); plot(val,fval(1:length(val))); xlabel('x'); ylabel('y'); fprintf(' 在 区 间 [%f,%f] 内 \n',x(k),x(k+1)); fprintf('三次样条函数 s（%d）=',k); pretty(ans) end

ቤተ መጻሕፍቲ ባይዱ
2, 3, ；
给出利用上述两种递推公式生成的序列的第 50 项。
解答：
第一小题： 程序： clc a(1) = sym(1); for j=1:49
a(j+1)=a(j)/sym(3); end %format rat disp('第五十项为：') disp(a(50)) 运行结果：
第二小题： 程序： clc a(1)=sym(1); a(2)=sym(1/3); for j=1:50
if i~=j w=w*(x0(j)-x0(i));
end end f(k)=f(k)+y0(j)/w; end end syms x px=f(1); for i=2:n d=1; for j=2:i d=d*(x-x0(j-1)); end px=px+f(i)*d; end p=px;
调用程序： function cdd_five a=-1; b=1; n=[3 6 9 12] for i=1:4 subplot(2,2,i) x=[a:0.001:b]; plot(x,f(x)); text(0,f(0),'y=f(x)'); hold on clear y0 y0=getxy(a,b,n(i)); p=cdd_Newton(a,b,y0,n(i)); fprintf('在 n=%d 时 Newton 插值多项式为： ',n(i)); px=simplify(p); px=vpa(px); pretty(px); plot(x,subs(p,x),'r'); text(0,subs(p,0),'y=p(x)'); end
pp(k)=simplify(pp(k)); coeff=sym2poly(pp(k)); ans=sym(coeff,'d'); ans=poly2sym(ans,'t'); val=x(k):0.001:x(k+1); for i=1:length(val)
fval(i)=coeff(1)*val(i)^3+coeff(2)*val(i)^2+c oeff(3)*val(i)+coeff(4);
function y=f(x) y=x^3+4*x^2-10; 运行结果：
Aitken 加速后： 程序：
clc; k=1; x(1)=1; while(f(x(k))~=0)
k=k+1; a=sqrt(10/(x(k-1)+4));
b=sqrt(10/(a+4)); x(k)=b-(b-a)^2/(b-2*a+x(k-1)); end fprintf('方程在[1,1.5]内的根为: x=%f 迭代次数为：k=%d',x(k),k); function varout=f(x) varout=x^3+4*x^2-10; 运行结果：
0.45 0.6708
0.53 0.7280
求解，其中边界条件为
.
三次样条差值函数程序： function cdd_chazhi(x,y,,s) n=size(x,2); for i=1:n-1
h(i)=x(i+1)-x(i); end for i=2:n-1
lm(i)=h(i)/(h(i)+h(i-1)); u(i)=1-lm(i);
z=z-A(n-i,j)*x(j); end x(n-i)=z/A(n-i,n-i); end fprintf('方程组的解为：\n'); x 运行结果：
题目：
四、已知一组数据点 数 满足
第四题
，编写一程序求解三次样条插值函
并针对下面一组具体实验数据
0.25
0.3
0.5000
0.5477
0.39 0.6245
程序：
function p=cdd_Newton(a,b,y0,n) for i=1:n
x0(i)=a+(i-1)*(b-a)/(n-1); %y0(i)=fx(x0(i)); %y0(i)=1/(1+25*x0(i)^2); end for k=1:n f(k)=0; for j=1:k
w=1; for i=1:k
g(i)=3*(u(i)*(y(i+1)-y(i))/h(i)+lm(i)*(y(i)-y(i1))/y(i-1)); end g(1)=3*(y(2)-y(1))/h(1)-h(1)*d2s(1)/2; g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)*d2s(2)/2 lm(n)=1; u(1)=1; lm u A=diag(lm(2:n),-1)+diag(u,1); for i=1:n
题目：
[0.39,0.45]
第五题
[0.45,0.53]
五、编写程序构造区间
上的以等分结点为插值结点的 Newton 插值公式，假设结点数
为 （包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出
相应的插值公式。以
，
为例计算其对应的插值公式，
分别取不同的 值并画出原函数的图像以及插值函数的图像，观察当 增大时的逼近效果.
function y=f(x) y=1./(1+25*x.^2); function y=getxy(a,b,n) for i=1:n
x0(i)=a+(i-1)*(b-a)/(n-1); y0(i)=f(x0(i)); end y=y0;
运行结果：
图形： 由结果可知，n 越大，插值函数曲线越逼近原函数
a(j+2)=10*a(j+1)/sym(3)-a(j); end format rat disp('第五十项为:') disp(a(50)) 运行结果：
第二题
题目：
二、利用迭代格式
xk 1
10 , k 0,1, 2, xk 4
及 Aitken 加速后的新迭代格式求方程 x3 4x2 10 0 在[1, 1.5] 内的根
迭代法计算停止的条件为：
max
1 j3
x (k 1) j

x(k) j
106 ．
2. 用 Gauss 列主元消去法、QR 方法求解如下方程组：
2

4
4
2 1 2
1 3 0
2 1 1


x1 x2 x3



1 2 1
QR 方法： 程序：
function three_QR(A) clc b=size(A); n=b(1); for i=1:n-1
Q=eye(n); a=A(i:n,i); a(1)=a(1)-sqrt(sum(a.^2)); w=a; H=eye(n-i+1)-2/(w'*w)*(w*w'); Q(i:n,i:n)=H; A=Q*A; end x=zeros(n,1); x(n)=A(n,n+1)/A(n,n); for i=1:n-1 z=A(n-i,n+1); for j=n-i+1:n;