《应用计算方法教程》matlab作业二.docx
- 文档编号:4422837
- 上传时间:2022-12-01
- 格式:DOCX
- 页数:15
- 大小:20.24KB
《应用计算方法教程》matlab作业二.docx
《《应用计算方法教程》matlab作业二.docx》由会员分享,可在线阅读,更多相关《《应用计算方法教程》matlab作业二.docx(15页珍藏版)》请在冰豆网上搜索。
《应用计算方法教程》matlab作业二
作业六
6-1试验目的计算特色值,实现算法
试验内容:
随机产生一个10阶整数矩阵,各数均在-5和5之间。
(1〕用MATLAB函数“eig〞求矩阵全部特色值。
(2〕用幂法求A的主特色值及对应的特色向量。
〔3〕
用根本QR算法求全部特色值〔可用
MATLAB函数“qr〞实现矩阵的
QR分解〕。
原理
幂法:
设矩阵A的特色值为|1|>|
2|
|
n|并设A有完好的特色向量系
1,
2,,n(它们线
性没关),那么对任意一个非零向量
V0
Rn所构造的向量序列
VkAVk1有lim
(Vk)j
1,其中
k
(Vk
1)j
(Vk)j表示向量的第
j个重量。
为防范逐次迭代向量
Vk不为零的重量变得很大〔|
1|
1
时〕或很小〔|
1|
1
时〕,将每一
步的Vk按其模最大的元素进行归一化。
详尽过程以下:
选择初始向量V0
,令mk
max(Vk),Uk
Vk
Vk1AUk,k1,当k
充分大时
mk
Uk
1
max(Vk1)
1。
max
(1)
QR法求全部特色值:
A
A1
Q1R1
RQ1
A2
Q2R2
1
k
1,2,3,
RkQk
Ak1
Qk1Rk1
由于此题的矩阵是10阶的,上述算法计算时间过长,考虑采用改良算法——移位加速。
迭代格式以下:
Ak
qkIQkRk
Ak
1RkQk
qkI
an(k
1,n)
1
an(k1,)
n
的特色值
(k)
(k)
,当
(k)
(k)
为实数时,选
qk为
计算Ak右下角的二阶矩阵
(k)
a
(k)
n1
n
n1
n
a
1
n,n
n,n
n(k
1),n(k)中最凑近an(k,n)的。
程序
A=-5+round(10*rand(10));
[V,D]=eig(A)
[lamdau]=lab6_2_power(A,[1;1;1;1;1;1;1;1;1;1],10^(-5),1000)
d=lab6_3_qr2(A,10^(-5))
function
lamda=0;
err=1;
[lamdau]=lab6_2_power(a,v,eps,N)
k=1;
while(k<=N&&err>eps)
u=a*v;
[mj]=max(abs(u));
dc=abs(lamda-m);
u=u/m;
dv=norm(u-v);
err=max(dc,dv);
v=u;
lamda=m;
k=k+1;
end
functionD=lab6_3_qr2(A,eps)
[n,n]=size(A);
m=n;
D=zeros(n,1);
B=A;
while(m>1)
while(abs(B(m,m-1))>=eps*(abs(B(m-1,m-1))+abs(B(m,m))))
S=eig(B(m-1:
m,m-1:
m));
[j,k]=min([abs(B(m,m)-S
(1)),abs(B(m,m)-S
(2))]);
[Q,U]=qr(B-S(k)*eye(m));
B=U*Q+S(k)*eye(m);
end
A(1:
m,1:
m)=B;
m=m-1;
B=A(1:
m,1:
m);
end
D=diag(A);
界面
(1)
(2)
(3)
作业七
7-1试验目的:
熟悉代数插值
试验内容:
在f(x)在7个点的函数值以下表所示,分别使用拉格朗日插值法和牛顿插值法求f(0.596)与f(0.906)的近似值。
xi
yi
1
原理
拉格朗日插值多项式:
n
(xx0)(xx1)(xxj1)(xxj1)
(xxn)
Ln(x)
j
0
(xj
x0)(xjx1)(xjxj1)(xjxj1)
(xjxn)yj
n
(
n
x
xi)yj
j
0
i
0
xj
xi
i
j
牛顿插值多项式:
Nn(x)
f[x0,
f(x0),xn](x
f[x0,x1](x
x0)(x
x0xn1
)
)
f[x0,x1,x2
](x
x0)(x
x1)
k
f(xj
)
其中f[x0,,xk]
。
j0(xjx0)
(xjxj1)(xjxj1)(xjxk)
程序
functiony1=lab7_1_Lagrange(x,y,x1)
y1=0;
[mn]=size(x);
n=n-1;
fork=1:
n+1
t=1;
fori=1:
n+1
if(i~=k)
t=t*(x1-x(i))/(x(k)-x(i));
end
end
y1=y1+t*y(k);
end
functiony1=lab7_2_Newton(x,y,x1)
[mn]=size(x);
n=n-1;
forj=1:
n
fori=n+1:
-1:
j+1
y(i)=(y(i)-y(i-1))/(x(i)-x(i-j));
end
end
y1=y(n+1);
forj=n:
-1:
1
y1=y(j)+(x1-x(j))*y1;
end
界面
作业八
8-1试验目的:
熟悉最小二乘法拟合多项式
试验内容:
给定数据点(xi,yi),
X
f(x)
用3次最小二乘多项式拟合数据,并求平方误差。
原理
要作三次最小二乘拟合,令
S(x)
a0
0(x)
a11(x)a22(x)
a3
3(x),
j(x)xj,计算法方程
(0,0)
(0,1)
(0,n)
a0
d1
Gad,其中G
(
1,0)
(
1,
1)
(
1,n),a
a1
,d
d2
,
(n,0)
(n,1)
(n,n)
an
dn
m
(
j,k)
(xi
)
j(xi)
k(xi),(j,k
0,1,,n)
i
0
(xi)
1。
m
dk
(f,k)
i
0
(xi)f(xi)k(xi),(k
0,1,
n)
其平方误差为s
m
f(xi)]2。
i0
(xi)[S(xi
)
程序
x=[0.40.550.650.800.901.05];
f=[0.410750.578150.696750.888111.026521.25386];
G=zeros(4,4);
forj=1:
4
fork=1:
4
fori=1:
6
G(j,k)=G(j,k)+x(i)^(j+k-2);
end
end
end
d=zeros(4,1);
fork=1:
4
fori=1:
6
d(k,1)=d(k,1)+f(i)*x(i)^(k-1);
end
end
a=G\d
s=0;
fori=1:
6
s=s+(f(i)-(a
(1)+a
(2)*x(i)+a(3)*x(i)^2+a(4)*x(i)^3))^2;
end
s
界面
作业九
9-1试验目的:
熟悉数值积分公式,掌握数值计算定积分的方法
试验内容:
采用不同样方法数值计算积分
1ln(1x)
dx
0
x
编写复合梯形公式和复合
Simpson公式通用子程序,分别采用
4,8,16,32,64均分区间
计算。
原理
复合梯形公式:
将区间
[a,b]作n均分,h
ba,结点xi
a
ih(0in),
n
Tn
n1h
h
n1
f(xi)
f(b)]
[f(xi)f(xi1)]
[f(a)
2
i02
2
i1
复合Simpson公式:
将区间[a,b]作2n均分,记h
b
a,
2n
S2n
n
1
h[f(x2i)4f(x2i1)f(x2i2)]
h[f(a)
4
n
1
f(x2i1)
2
n
1
f(x2i)f(b)]
i
0
3
3
i
0
i
1
程序
functiony=lab9_f(x)
y=(log(1+x))/x;
functiony=lab9_1_fTrapezoid(a,b,eps,n)
f0=0;
h=(b-a)/n;
fori=0:
n-1
f0=f0+lab9_f(a+i*(b-a)/n)+lab9_f(a+(i+1)*(b-a)/n);
end
y=f0*h/2;
functiony=lab9_2_fSimpson(a,b,eps,n)
f0=0;
h=(b-a)/n;
k=n/2;
fori=0:
k-1
f0=f0+lab9_f(a+2*i*(b-a)/n)+4*lab9_f(a+(2*i+1)*(b-a)/n)+lab9_f(a+(2*i+2)*(b-a)/n);
end
y=f0*h/3;
界面
作业十
10-1试验目的:
学会用Euler法、改良Euler法、经典的4阶Runge-Kutta法求解常微分方程
初值问题。
试验内容:
分别用
1)Euler法〔步长〕
2)改良Euler法〔步长〕
3)4阶Runge-Kutta〔步长〕
求解下面的初值问题:
y
(y
1)(y
3),0x
2
比较公共节点解的误差。
精确解为
y
3
2(1
e2x)1。
原理
Euler法:
令f(x,y(x))
y(x),yn
1
yn
hf(xn
yn)〔n
ba〕。
h
k1
f(xn,yn)
改良Euler
法〔梯形公式〕:
k2
f(xn
h,yn
hk1)
y
1
y
n
h(k
k
2
)
n
1
2
k1
hf(xn,yn)
k2
hf(xn
h
yn
1
k1)
2
2
4阶Runge-Kutta:
k3
hf(xn
h,yn
1k2)
2
2
k4
hf(xn
h,yn
k3)
yn1
yn
1
2k2
2k3
k4)
(k1
6
程序
y0=-2;
h=0.025;
n=2/h;
y
(1)=y0-h*(y0+1)*(y0+3);
fori=2:
n
y(i)=y(i-1)-h*(y(i-1)+1)*(y(i-1)+3);
end
y
forj=1:
n
y1(j)=-3+2/(1+exp(-4*j/n));
end
y1
x=0.025:
0.025:
2;
plot(x,y1,'r'
x,y,
'b'
)
y0=-2;
h=0.05;
n=2/h;
k1=-(y0+1)*(y0+3);
k2=-(y0+h*k1+1)*(y0+h*k1+3);
y21
(1)=y0+h/2*(k1+k2);
fori=2:
n
k1=-(y21(i-1)+1)*(y21(i-1)+3);
k2=-(y21(i-1)+h*k1+1)*(y21(i-1)+h*k1+3);
y21(i)=y21(i-1)+h/2*(k1+k2);
end
y21
forj=1:
n
y22(j)=-3+2/(1+exp(-4*j/n));
end
y22
x=0.05:
0.05:
2;
plot(x,y21,
'r'
x,y22,
'b'
)
y0=-2;
h=0.1;
n=2/h;
k1=-h*(y0+1)*(y0+3);
k2=-h*(y0+k1/2+1)*(y0+k1/2+3);
k3=-h*(y0+k2/2+1)*(y0+k2/2+3);
k4=-h*(y0+k3+1)*(y0+k3+3);
y31
(1)=y0+(k1+2*k2+2*k3+k4)/6;
fori=2:
n
k1=-h*(y31(i-1)+1)*(y31(i-1)+3);
k2=-h*(y31(i-1)+k1/2+1)*(y31(i-1)+k1/2+3);
k3=-h*(y31(i-1)+k2/2+1)*(y31(i-1)+k2/2+3);
k4=-h*(y31(i-1)+k3+1)*(y31(i-1)+k3+3);
y31(i)=y31(i-1)+(k1+2*k2+2*k3+k4)/6;
end
y31
forj=1:
n
y32(j)=-3+2/(1+exp(-4*j/n));
end
y32
x=0.1:
0.1:
2;
plot(x,y31,
'r'
x,y32,
'b'
)
界面
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 应用计算方法教程 应用 计算方法 教程 matlab 作业
![提示](https://static.bdocx.com/images/bang_tan.gif)