常州大学数值分析课后习题答案第二章第三章第四章节资料.docx
- 文档编号:12149231
- 上传时间:2023-04-17
- 格式:DOCX
- 页数:16
- 大小:29.78KB
常州大学数值分析课后习题答案第二章第三章第四章节资料.docx
《常州大学数值分析课后习题答案第二章第三章第四章节资料.docx》由会员分享,可在线阅读,更多相关《常州大学数值分析课后习题答案第二章第三章第四章节资料.docx(16页珍藏版)》请在冰豆网上搜索。
常州大学数值分析课后习题答案第二章第三章第四章节资料
数值分析作业
第二章
1、用Gauss消元法求解下列方程组:
2x1-x2+3x3=1,
(1)4x1+2x2+5x3=4,
x1+2x2=7;
(2)解:
A=[2-131;4254;1207]
n=size(A,1);x=zeros(n,1);flag=1;
%消元过程
fork=1:
n-1
fori=k+1:
n
ifabs(A(k,k))>eps
A(i,k+1:
n+1)=A(i,k+1:
n+1)-A(k,k+1:
n+1)*A(i,k)/A(k,k);
else
flag=0;
return
end
end
end
%回代过程
ifabs(A(n,n))>eps
x(n)=A(n,n+1)/A(n,n);
else
flag=0;
return
end
fori=n-1:
-1:
1
x(i)=(A(i,n+1)-A(i,i+1:
n)*x(i+1:
n))/A(i,i);
end
return
x
A=2-131
4254
1207
x=9
-1
-6
11x1-3x2-2x3=3,
(2)-23x1+11x2+1x3=0,
x1+2x2+2x3=-1;
(2)解:
A=[11-3-23;-231110;122-1]
n=size(A,1);x=zeros(n,1);flag=1;
%消元过程
fork=1:
n-1
fori=k+1:
n
ifabs(A(k,k))>eps
A(i,k+1:
n+1)=
A(i,k+1:
n+1)-A(k,k+1:
n+1)*A(i,k)/A(k,k);
else
flag=0;
return
end
end
end
%回代过程
ifabs(A(n,n))>eps
x(n)=A(n,n+1)/A(n,n);
else
flag=0;
return
end
fori=n-1:
-1:
1
x(i)=(A(i,n+1)-A(i,i+1:
n)*x(i+1:
n))/A(i,i);
end
return
x
A=11-3-23
-231110
122-1
x=0.2124
0.5492
-1.1554
4、用Cholesky分解法解方程组
323x15
220x23
3012x37
解:
.
A=[323;220;3012];
b=[537];
lambda=eig(A);
iflambda>eps&isequal(A,A')
[n,n]=size(A);
R=chol(A);
%解R'y=b
y
(1)=b
(1)/R(1,1);
ifn>1
fori=2:
n
y(i)=(b(i)-R(1:
i-1,i)'*y(1:
i-1)')/R(i,i);
end
end
%解Rx=y
x(n)=y(n)/R(n,n);
ifn>1
fori=n-1:
-1:
1
x(i)=(y(i)-R(i,i+1:
n)*x(i+1:
n)')/R(i,i);
end
end
x=x';
else
x=[];
disp('该方法只适用于对称正定的系数矩阵!
');
end
R=1.73211.15471.7321
00.8165-2.4495
001.7321
y=2.8868-0.40820.5774
x=1.00000.50000.3333
5.用列主元Doolittle分解法解方程组
解:
A=[345;-134;-23-5;];345X12
b=[2,-26]';-134X2-2
[L,U,pv]=luex(A);-23-5X36
y=L\b(pv);
x=U\y
结果如下:
x= 1
1
-1
14.已知,计算
.
解:
A=[10099;9998];
cond(A,inf)
ans=3.9601e+04
cond(A,2)
ans=3.9206e+04
27.编写LU分解法,改进平方根法,追赶法的Matlab程序,并进行相关数值试验。
解:
LU分解法程序
Function[L,U]=lup(A)
%lup:
LUfactorization
%Synopsis:
[L,U]=lup(A)
%Input:
A=coefficientmatrix
%Output:
L:
lowertriangularmatrix
%Uuppertriangularmatrix
Formatshort
[m,n]=size(A);
Ifm~=n,error(`Amatrixneedstobesquare`);
End
Pv=(1:
n)`;
%LUfactorization
Fori=1:
n-1
Pivot=A(i,i);
Fork=i+1;n
A(k,i)=A(k,i)/pivot;
A(k,i+1;n)=A(k,i+1;n)-A(k,i)*A(i,i+1;n);
End
End
L=eye(size(A))+tril(A,-1);
%extractLandU
U=triu(A)
改进平方根法程序
Function[x]=ave(A,b,n)
L=zeros(n,n);
D=diag(n,0);
S=L*D;
Fori=1:
n
L(i;i)=1;
End
Fori=1;n
Forj=1;n
If(eig(A)<=0)|(A(i,j)~=A(j,i))disp(`wrong`);
Break;
End
End
End
D(1,1)=A(1,1);
Fori=2;n
Forj=1;i-1
S(i,j)=A(i,j)-sum(S(i,1;i-1)*L(j,1;j-1)`);
L(i,1;i-1)=S(i,1;i-1)/D(1;i-1,1;i-1);
end
D(i,i)=A(i,i)-sum(L(i,1;i-1)*L(i,1;i-1)`);
End
D(i,i)=A(i,i)-sum(S(i,1;i-1)*D(1;i-1,1;i-1)*y(1;i-1)))/D(i,i);
End
Y=zeros(n,1);
X=zero(n,1);
Fori=1;n
Y(i)=(b(i)-sum(L(i,1;i-1)*D(1;i-1,1;i-1)*y(1;i-1)))/D(i,i);
End
Fori=n;-1;1
X(i)=y(i)-sum(L(i+1;n,i)`*x(i+1;n));
End
追赶法程序
Function[x,L,U]=Thomas(a,b,c,f)
N=length(b);
%对A进行分解
U
(1)=b
(1);
Fori=2;n
If(u(i-1)`=0)
L(i-1)=a(i-1)/u(i-1);
U(i)=b(i)-l(i-1)*c(i-1);
Else
Break;
End
End
L=eye(n)+diag(1,-1);
U=diag(u)+diag(c,1);
X=zeros(n,1);
Y=x;
%?
求解ly=b?
Y
(1)=f
(1);
Fori=2;n
Y(i)=f(i)-l(i-1)*y(i-1);
End
%?
求解Ux=y?
If(u(n)`=0)
X(n)=y(n)/u(n);
End
Fori=n-1;-1;1
X(i)=(y(i)-c(i)*x(i+1))/u(i);
End
第三章
1、设节点x0=0,x1=π/8,x2=π/4,x3=3π/8,x4=π/2,适当选取上述节点用Lagrange插值法分别构造cosx在区间[0,π/2]上的一次,二次和四次插值多项式P1(x)P2(x)和P4(x),并分别计算P1(x),P2(x),P4(x)其中X取π/3。
A=fliplr(A);Return
x=[π/8,3π/8];y=cos(x);x0=π/3;
[A,Y]=lagrange(x,y,x0);P1=vpa(poly2sym(A),3)Y
P1=1.19x-0.689Y=0.4729x0=π/3;
[A,Y]=lagrange(x,y,x0);P2=vpa(poly2sym(A),3)Y
P2=x2-0.109x-0.336Y=0.5174
x=[0,π/8,π/4,3π/8,π/2];y=cos(x);x0=π/3;
[A,Y]=lagrange(x,y,x0);P4=vpa(poly2sym(A),3)Y
P4=x4+0.00282x3-0.514x2+0.0232x+0.0287Y=0.5001
7.根据列表函数
x
0.0
1.0
2.0
3.0
4.0
f(x)
0.50
1.25
2.75
3.50
2.75
选取适当的节点,用逐次线性插值法给出三次多项式在2.8处的值。
答:
Matlab程序function
[T,y0]=aitken(x,y,x0,T0)ifnargin==3T0=[];
end
n0=size(T0,1);
m=max(size(x));
n=n0+m;
T=zeros(n,n+1);
T(1:
n0,1:
n0+1)=T0;
T(n0+1:
n,1)=x;
T(n0+1:
n,2)=y;
ifn0==0i0=2;
else
i0=n0+1;
end
fori=i0:
n
forj=3:
i+1
T(i,j)=fun(T(j-2,1),T(i,1),T(j-2,j-1),T(i,j-1),x0);
end
y0=T(n,n+1);
return
function[y]=fun(x1,x2,y1,y2,x)y=y1+(y2-y1)*(x-x1)/(x2-x1);
return
%选取0、1、3、4四个节点,求三次插值多项式x=[0,1,3,4];
y=[0.5,1.25,3.5,2.75];
x0=2.8;
[T,y0]=aitken(x,y,x0)
y0=3.419000000000000
8.根据上题中的列表函数,写出差商表,并写出Newton插值多项式N2(x)和N4(x)。
答:
差商表:
x
f(x)
0.0
0.5
1.0
1.25
0.75
2.0
2.75
1.125
0.375
3.0
3.50
1
0.125
-0.25
4.0
2.75
0.5625
-0.0625
-0.21875
0.03125
由Nn(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)+…+an(x-x0)(x-x1)…(x-xn-1)得
N2(x)=0.5+0.75(x-0.00)+0.375(x-0.0)(x-1.0)
=0.375x2+0.375x+0.5
N4(x)=0.5+0.75(x-0.00)+0.375(x-0.0)(x-1.0)-0.25(x-0.0)(x-1.0)(x-2.0)+0.03125(x-0.0)(x-2.0)(x-1.0)(x-3.0)
=0.03125x4-0.4375x3+1.46875x2-0.3125x+0.5
16、选取适当的函数y=f(x)和插值节点,编写Matlab程序,分别利用Lagrange插值方法,Newton插值方法确定的插值多项式,并将函数y=f(x)的插值多项式和插值余项的图形画在同一坐标系中,观测节点变化对插值余项的影响。
答:
Matlab程序function
[C,D,Y]=newpoly(x0,y0,x)
%检验输入参数
ifnargin<2|nargin>3error('IncorrectNumberofInputs');end
iflength(x0)~=length(y0)error('Thelengthofx0mustbeequaltoitofy0');end
n=length(x0);D=zeros(n,n);D(:
1)=y0';
%计算差商表
forj=2:
nfork=j:
n
ifabs(x0(k)-x0(k-j+1)) D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k)-x0(k-j+1));endend %计算Newton插值多项式的系数C=D(n,n); fork=(n-1): -1: 1 C=conv(C,poly(x0(k)));m=length(C); C(m)=C(m)+D(k,k);end ifnargin==3 Y=polyval(C,x);end x=[01234]; y=[0.5,1.25,2.75,3.5,2.75];x0=[01234]; y0=[0.5,1.25,2.75,3.5,2.75];%用lagrang插值法计算[A,Y]=lagrange(x,y,x0) Lx=vpa(poly2sym(A),4)%用newton插值法计算[C,D,X]=newpoly(x0,y0,x)Nx=vpa(poly2sym(C),4)%绘制两者图像plot(x,Y,'b*',x0,X,'r-') A=0.5000-0.31251.4687-0.43750.0313 Y=0.50001.25002.75003.50002.7500 Lx=0.5x4-0.3125x3+1.469x2-0.4375x+0.03125 C=0.5000-0.31251.4688-0.43750.0313 D=0.50000000 1.25000.7500000 2.75001.50000.375000 3.50000.7500-0.3750-0.25000 2.7500-0.7500-0.7500-0.12500.0313 X=0.50001.25002.75003.50002.7500 Nx=0.5x4-0.3125x3+1.469x2-0.4375x+0.0312 第四章 6、已知列表函数 x1.02.03.04.05.0 y1.2222.9845.4668.90213.592 用最小二乘法求形如y=axebx的拟合函数。 答: Matlab程序 function[a,b]=ec(x,y)Y=log(y)'; A=zeros(5,3); fori=1: 5A(i,1)=1; A(i,2)=log(x(i)); A(i,3)=i;end c=inv(A'*A)*(A'*Y); a=exp(c (1)); b=c(3); fori=1: 5 y=a*x.*exp(b*x); endreturn x=[12345]; y=[1.2222.9845.4668.90213.592]; [a,b]=ec(x,y)输出结果为: a=1.000202219673205b=0.200293860504786 8、学习Matlab内部的函数lsqcurvefit,并设计数值实验使用lsqcurvefit。 答: Matlab内部函数lsqcurvefit是用来解决非线性拟合的最小二乘问题的。 其调用格式为: x=lsqcurvefit(fun,x0,xdata,ydata)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub) x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options) [x,resnorm]=lsqcurvefit(…) [x,resnorm,residual]=lsqcurvefit(…) [x,resnorm,residual,exitflag]=lsqcurvefit(…) [x,resnorm,residual,exitflag,output]=lsqcurvefit(…) [x,resnorm,residual,exitflag,output,lambda]=lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(…) 输入参数: fun为待拟合函数,计算x处拟合函数值,其定义为functionF=myfun(x,xdata) x0为初始解向量,即拟合参数的初始解;xdata,ydata为满足关系ydata=F(x,xdata)的数据; lb、ub为解向量的下界和上界lb≤x≤ub,若没有指定界,则lb=[],ub=[];options为指定的优化参数;输出参数: x为迭代得出解向量,即拟合出的参数;resnorm=sum((fun(x,xdata)-ydata).^2),即x处残差平方和,最小二乘式值; residual=fun(x,xdata)-ydata,即在x处的残差; exitflag为终止迭代的条件;output为输出的优化信息; lambda为解x处的Lagrange乘子; jacobian为解x处拟合函数fun的jacobian矩阵。 functionF=myfun(x,xdata) F=(x (1).*xdata).*(exp(x (2).*xdata)); xdata=[1,2,3,4,5]; ydata=[1.222,2.984,5.466,8.902,13.592];x0=[0,0]; [x,resnorm]=lsqcurvefit(@myfun,x0,xdata,ydata) 输出结果为: Localminimumfound. Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueofthefunctiontolerance. 0.9999583489763910.200014132812834resnorm=8.067930437509675e-7
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 常州 大学 数值 分析 课后 习题 答案 第二 第三 第四 章节 资料