数值分析上机题目

s=s+f(a+(i-1/2)*h);
end
T(col,1)=1/2*T(col-1,1)+h/2*s;
%以上为第一列，也就是半分区间得到的复化梯形公式求积
else
T(col+1-j,j)=(4^(j-1)*T(col+2-j,j-1)-T(col+1-j,j-1))/(4^(j-1)-1);
%以上为第2到col列，也就是第1到第col-1次外推公式
a.Matlab代码，five.m
%矩阵的三角分解 P142 练习五第6题
clear();
format long;
a=[2 3 4 0;1 1 9 2;1 2 -6 1];
for j=1:length(a)-1
c=max(abs(a(j:length(a)-1,j)));
[x,y]=find(abs(a(j:length(a)-1,j)==c));
y(:,1)=[1.25,2.50,1.00,5.50]';%已知点集合x和y
symst w;
w(1)=1;
%计算基函数序列w和差商表y，以及函数序列的权数diag(y),计算的牛顿三次多项式表述为t的函数
for j=2:length(x)
fori=j:length(x)
y(i,j)=(y(i,j-1)-y(i-1,j-1))/(x(i)-x(i-j+1));
%%%%%
x(:,1)=[0,0,0]';
for col=2:3
for row=1:length(A)
if row来自百度文库=1
x(row,col)=1/D(row,row)*(b(row)+U(row,row+1:length(A))*x(row+1:length(A),col-1));
elseifrow==length(A)
i=i+1;
end
w(j)=prod(t-x(1:j-1));
j=j+1;
end
disp('三次牛顿插值多项式为');
disp(collect(w*diag(y)));
plot(x,y(:,1),'*');
hold on;
fplot(collect(w*diag(y)),[-0.5,2.5]);
legend({'已知点集','三次牛顿插值多项式函数'},'location','northwest','fontsize',14);
end
end
end
disp('从第2列开始是外推得到的结果，第1列是使用复化梯形公式得到的结果');
disp(T);
b.Matlab函数文件，f.m
function y=f(x)
if x==0
y=1;
else
y=sin(x)/x;
end
c.运行结果如下：
四．
L5.6用Gauss-Jordan消去法解方程组
end
end
end
while norm(X(:,col)-X(:,col-1),inf)>10^-5
col=col+1;
for row=1:length(A)
if row==1
X(row,col)=1/D(row,row)*(b(row)+U(row,row+1:length(A))*X(row+1:length(A),col-1));
clear();
format long;
A=[1 1 0.5;1 1 0.25;0.5 0.25 2;];
r0=2.5365258;%主特征值的准确值（保留八位有效数字）
v(:,1)=[1 1 1]';%v向量中包含主特征值的序列，其最大值收敛到主特征值
u(:,1)=[1 1 1]';%u为主特征向量的序列
Gn(i,j)=int(f(i)*f(j),-1,1);
j=j+1;
end
b(i)=int(f(i)*y,-1,1);
i=i+1;
end
a=b/Gn;%最佳平方逼近的系数行向量
disp('逼近函数表达式');
disp(vpa(f*a'));
disp('最佳函数逼近得平方误差');
disp(vpa(int(y^2,-1,1)-a*b'));
if x~=1
d=a(j,:);
a(j,:)=a(x+j-1,:);
a(x+j-1,:)=d;
end
l=a(:,j)/a(j,j);
%l=l(j+1:end);
fori=1:length(a)-1
ifi~=j
a(i,:)=a(i,:)-a(j,:)*l(i);
end
i=i+1;
end
end
disp('经过转化后的系数矩阵');
数值分析上机编程题目
1631110XXXX
材料科学与工程学院
一
L2.7给定数据表2-15.用Newton插值公式计算3次插值多项式N3(x).
表2-15
x
1
1.5
0
2
f(x)
1.25
2.50
1.00
5.50
a.Matlab代码如下，two.m,
%第二章，P45,练习题2第七题
clear();
x=[1,1.5,0,2];
elseifrow==length(A)
X(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*X(1:row-1,col));
else
X(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*X(1:row-1,col)+U(row,row+1:length(A))*X(row+1:length(A),col-1));
end
end
end
disp('Gauss-Seidel迭代使得结果向量的无穷范数不超过10^-5，第一行为初始值[0,0,0]');
disp(X');
b.运行结果如下：
六．
例7.9 用Steffensen迭代法求方程f(x)=x3-x-1=0的实根。
a.Matlab代码，seven.m
%第7章，非线性方程的迭代解法，P180，例7.9，Steffensen迭代法
end
end
end
while norm(x(:,col)-x(:,col-1),inf)>10^-5
col=col+1;
for row=1:length(A)
if row==1
x(row,col)=1/D(row,row)*(b(row)+U(row,row+1:length(A))*x(row+1:length(A),col-1));
clear();
symsx;
%所使用的非线性基函数序列，用符号表示
y=abs(x);%被逼近函数
f=[1,x^2,x^4];
%求解法方程的系数矩阵a*Gn=b，其中a和b均为行向量
Gn=ones(length(f),length(f));
fori=1:length(f)
for j=1:length(f)
%T(1,col)对应的计算的是多少步的值,col→coln关系
col=col+1;%此时求得是第n+1次均分后的结果，使用的是第n次的结果，注意在矩阵
%计算的第n斜列是第n-1次均分的结果
for j=1:col
if j==1
h=(b-a)/2^(col-2);%使用n+1之前的第n次结果
s=0;
fori=1:2^(col-2)
L=zeros(length(A));
U=zeros(length(A));
for col=1:length(A)-1
for row=col+1:length(A)
L(row,col)=-A(row,col);
row=row+1;
end
col=col+1;
end
U=-A-L+D;
%
%Jacobi迭代法的使用！！！
end
end
end
disp('Jocobi迭代使得结果向量的无穷范数不超过10^-5，第一行为初始值[0,0,0]');
disp(x');
%
%Gauss-Seidel迭代法的使用！！！
%
X(:,1)=[0,0,0]';
for col=2:3
for row=1:length(A)
if row==1
X(row,col)=1/D(row,row)*(b(row)+U(row,row+1:length(A))*X(row+1:length(A),col-1));
elseifrow==length(A)
X(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*X(1:row-1,col));
else
X(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*X(1:row-1,col)+U(row,row+1:length(A))*X(row+1:length(A),col-1));
ylabel('y','FontSize',20,'FontWeight','bold');
hold off;
b.运行结果如下：
三.
例4.9用Romberg求积法计算定积分 ，使计算值的误差不超过ε<0.5 10-6。
a.Matlab代码，four.m
%Romberg求积公式，外推原理
clear();
i=i+1;
x(i)=F(x(i-1));
end
disp('Steffensen迭代法，在误差小于10^-7时的迭代结果');
disp(x');
b.运行结果如下：
七．
例9.3 用幂法求矩阵
A=
的主特征值和主特征向量。
a.Matlab代码，nine.m
%第九章，P212，例9.3，使用幂法求主特征值和主特征向量
b.运行结果如下：
五．
例6.2 用J法和GS法分别求解方程组。
a.Matlab代码，six.m
%第六章,P149 例6.2 线性方程组的迭代解法
format long;
clear();
A=[10 -2 -2;-2 10 -1;-1 -2 3];
b=[1 0.5 1]';
D=diag(diag(A));
disp(a); %进行转化后的系数矩阵
%回代过程
X=zeros(1,length(a)-1);
fori=length(a)-1:-1:1
X(i)=(a(i,length(a))-X*a(i,1:length(a)-1)')/a(i,i);
i=i+1;
end
disp('求解的结果');
disp(X);
clear();
f=@(x)x^3-1;%方程f(x)=x^3-x-1=0的迭代公式可以选择为x=x^3-1,虽然是发散的
F=@(x)x-(f(x)-x)^2/(f(f(x))-2*f(x)+x);
x(1)=1.5;
x(2)=F(x(1));
i=2;
while abs(x(i)-x(i-1))>10^-7
fplot(y,[-1,1]);
hold on;
fplot(a*f',[-1,1]);
legend({'被逼近函数','逼近函数'},'Location','north','Orientation','horizontal','FontSize',16,'FontWeight','bold');
xlabel('x','FontSize',20,'FontWeight','bold');
clear();
format long;
a=0;
b=1;
T(1,1)=(b-a)/2*(f(a)+f(b));
T(2,1)=1/2*T(1,1)+(b-a)/2*f((a+b)/2);
T(1,2)=(4*T(2,1)-T(1,1))/(4-1);
col=2;
while abs(T(1,col)-T(1,col-1))>0.5*10^-6
xlabel('x','fontsize',16);
ylabel('y','fontsize',16);
hold off;
b.计算结果如下：
二．
L3.2 设f(x)=|x|，x∈[-1,1]。求在Φ=span{1,x2,x4}上的最佳平方逼近。
a.Matlab代码，three.m，
%第三章函数逼近与数据拟合，P68练习题，第2题
x(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*x(1:row-1,col-1));
else
x(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*x(1:row-1,col-1)+U(row,row+1:length(A))*x(row+1:length(A),col-1));
elseifrow==length(A)
x(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*x(1:row-1,col-1));
else
x(row,col)=1/D(row,row)*(b(row)+L(row,1:row-1)*x(1:row-1,col-1)+U(row,row+1:length(A))*x(row+1:length(A),col-1));