精品大连理工矩阵上机作业Word格式.docx
《精品大连理工矩阵上机作业Word格式.docx》由会员分享，可在线阅读，更多相关《精品大连理工矩阵上机作业Word格式.docx（18页珍藏版）》请在冰点文库上搜索。
w=[1,2,3,4,5,6,7,8,9,10];
d1=0;
d2=0;
d3=0;
y1=polyfit（m,w,1）;
y2=polyfit（m,w,2）;
y3=polyfit（m,w,3）;
y2=poly2sym（s2）;
y3=poly2sym（s3）;
y4=poly2sym（s4）;
f1=subs（y1,17）;
f2=subs（y2,17）;
f3=subs（y3,17）;
forp=1:
10;
d1=d1+（subs（y1,w（p））-m（p））^2;
d2=d2+（subs（y2,w（p））-m（p））^2;
d3=d3+（subs（y3,w（p））-m（p））^2;
d1=sqrt（d1）;
d2=sqrt（d2）;
d3=sqrt（d3）;
fun=[f1f2f3;
d2d3d4];
return;
结果
三次函数的均方误差最小，拟合的最好。
函数
functionf=fun（x）
symsa
a=x;
f=a*a*a+a*a+a-3;
梯度函数
functiondf=dfun（x）
df=3*x*x+2*x+1;
Newton法
functionresult=didainewton（x0）
k=0;
xk=x0;
xi=1;
e0=abs（x0-xi）;
ek=e0;
m=zeros（7,1）;
n=zeros（7,1）;
p=zeros（7,1）;
result=zeros（7,3）;
whilek<
7
ak=feval（'
fun'
xk）;
bk=feval（'
dfun'
xk=xk-ak/bk;
e0=ek;
k=k+1;
m（k）=xk;
ek=abs（m（k）-xi）;
jingdu=ek/（e0*e0）;
n（k）=ek;
p（k）=jingdu;
result=[m,n,p];
计算结果
i
Xi
Ei=|xi-x*|
Ei/Ei-12
-0.7
1.7
_
1
2.6206
1.6206
0.5607
2
1.7084
0.7084
0.2698
3
1.2064
0.2064
0.4112
4
1.0242
0.0242
0.5673
5
1.0004
0.0004
0.6535
6
1.0000
0.0000
0.6665
0.6609
Newton迭代法的收敛速度快，至少是平方收敛的，
GAUSS消去法
functionx=DelGauss（N）
%GaussÏ
û
È
¥
·
¨
symsM;
M=N;
a=zeros（M）;
b=ones（M,1）;
M
a（i,j）=1/（i+j-1）;
[n,m]=size（a）;
nb=length（b）;
det=1;
x=zeros（n,1）;
fork=1:
n-1
fori=k+1:
n
ifa（k,k）==0
return
m=a（i,k）/a（k,k）;
forj=k+1:
a（i,j）=a（i,j）-m*a（k,j）;
b（i）=b（i）-m*b（k）;
det=det*a（k,k）;
det=det*a（n,n）;
fork=n:
-1:
1
b（k）=b（k）-a（k,j）*x（j）;
x（k）=b（k）/a（k,k）
N=4
X=[-4.000060.0000-180.0000140.0000]T
N=8
X=[-8.0000001,504.00001,-7560.0,46200.0,-138600.0,216216.0,
-168168.0,51480.0]T
Jacobbi
functionx=jacobi（N）
symsM
A=zeros（M）;
A（i,j）=1/（i+j-1）;
p=zeros（M,1）;
q=zeros（M,1）;
q（i）=A（i,i）;
D=diag（q）;
L=-tril（A,-1）;
U=-triu（A,1）;
B=（inv（D））*（L+U）;
f=（inv（D））*b;
x=B*p+f;
n=1;
sigmal=1e-8;
whilenorm（（p-x））>
=sigmal
p=x;
x=B*p+f;
n=n+1;
eval（'
x'
N=4时
R（B）=2.5821>
1，迭代法不收敛
N=8时
R（B）=6.0422>
1,迭代法不收敛
Gauss-Seidel
functionx=gaussseidel（N）
D=diag（diag（A））;
G=（inv（D-L））*U;
f=（inv（D-L））*b;
x=G*p+f;
detal=1e-4;
whileabs（x-p）>
=detal
x=G*p+f;
eval（'
X=[2.8527-14.7133-3.538826.7350]T
X=[-3.920228.0256-9.4820-41.3197-34.0806-6.447528.846265.3427]T
共轭梯度法
functionresult=cg（N）
x0=zeros（M,1）;
r0=b-A*x0;
p0=r0;
sigmal=1e-4;
rk=r0;
pk=p0;
whiledot（r0,r0）>
=sigmal
ak=（dot（r0,r0））/（dot（pk,A*pk））;
ri=r0;
xk=x0+ak*pk;
rk=r0-ak*A*pk;
bk=（dot（rk,rk））/（dot（ri,ri））;
pk=rk+bk*p0;
p0=pk;
r0=rk;
x0=xk;
result=x0;
X=[-4.000060.0000-180.0000140.0000]T
X=[4.5773-72.4211199.8408-9.9954-174.8581-171.8511-10.8652269.4734]T
第五题
functionf=san（x1,y）
n=length（x1）;
h=zeros（1,n-1）;
n-1
h（p）=x1（p+1）-x1（p）;
n-2
r（p）=h（p+1）/（h（p+1）+h（p））;
u（p+1）=h（p）/（h（p+1）+h（p））;
u
（1）=1;
r（n-1）=1;
c=2*ones（1,n）;
A=diag（r,-1）+diag（u,1）+diag（c）;
g=zeros（n,1）;
g
（1）=3*（y
（2）-y
（1））/h
（1）;
g（n）=3*（y（n）-y（n-1））/h（n-1）;
forp=1;
n-2;
g（p+1）=3*（（y（p+2）-y（p+1））/h（p+1）*u（p+1）+r（p）*（y（p+1）-y（p））/h（p））;
m=A\g;
m=m'
f=m;
return
x1=[0:
4];
y=[1,3,3,4,2];
san（x1,y）
ans=
2.50001.0000-0.50001.0000-3.5000
利用教科书p179的公式（5-35）
functionfun=sanci（m,h,xk,y）
symsx
n=length（m）;
l=length（xk）;
fun=（（h（i）+2*（x-xk（i）））*y（i）/h（i）+（x-xk（i））*m（i））*（x-xk（i+1））^2/h（i）^2+（（h（i）-2*（x-xk（i+1）））*y（i+1）/h（i）+（x-xk（i+1））*m（i+1））*（x-xk（i））^2/h（i）^2
0<
=x<
（（9*x）/2+1）*（x-1）^2-x^2*（5*x-8）
1<
（7*x-4）*（x-2）^2-（（13*x）/2-16）*（x-1）^2
2<
（（11*x）/2-8）*（x-3）^2-（7*x-25）*（x-2）^2
3<
（9*x-23）*（x-4）^2-（（15*x）/2-32）*（x-3）^2
第六题
symsa;
f=a*a*a+2*a*a+10*a-100;
functiongrad=gfun（x）
grad=3*a*a+4*a+10;
迭代函数
functionresult=xianjiefa（x0）
sigmal=1e-6;
maxi=500;
m=zeros（500,1）;
d0=fun（x0）/gfun（x0）;
x1=x0-d0;
k=1;
xk=x1;
m
（1）=x1;
maxi
dk=fun（xk）*（xk-x0）/（fun（xk）-fun（x0））;
x0=xk;
xk=xk-dk;
ifabs（xk-x0）<
sigmal
break;
n=zeros（k,1）;
k
n（i）=m（i）;
result=（n）;
X0=0
X0=1
X0=2
X0=3
X0=4
X0=5
10.0000
6.1176
4.1333
3.5102
3.5135
3.8095
0.7692
2.2649
3.2399
3.4556
3.4661
3.5495
1.4177
2.9318
3.4316
3.4605
3.4606
3.4667
5.788
3.6259
3.4620
3.4607
2.5766
3.4417
3.1128
3.4600
3.5356
3.455
第七题
functionfun=intgra（n）
symsx;
symsm;
sum1=0;
sum2=0;
m=n;
f=（exp（3*x））*cos（pi*x）;
x0=0:
2*pi/m:
2*pi;
a=0;
b=2*pi;
h=（b-a）/m;
m-1
xk=a+i*h;
sum1=sum1+subs（f,xk）;
fori=0:
xk=a+（i+1/2）*h;
sum2=sum2+subs（f,xk）;
t=（（b-a）/（2*n））*（subs（f,a）+2*sum1+subs（f,b））;
s=（（b-a）/（6*n））*（subs（f,a）+2*sum1+4*sum2+subs（f,b））;
z=int（f,x,0,2*pi）;
T=vpa（t,2）;
S=vpa（s,2）;
Z=vpa（z,2）
Dt=abs（T-Z）;
Ds=abs（S-Z）;
fun=[T;
S;
Dt;
Ds];
第八题
Gauss-Chebyshev函数
functionresult=chebyshev（f,n）
x=zeros（1,n）;
A=zeros（1,n）;
m=n-1;
m
x（i+1）=cos（（2*i+1）*pi/（2*n））;
A（i+1）=pi/n;
result=double（sum（A（1:
n）.*subs（f,x（1:
n））））;
Gauss-Legendre函数
functionresult=lanendre（f,a,x）
t=pi/4+pi/4*x;
result=double（pi/4*sum（a.*subs（f,t）））;
计算过程
f1=x^2
f2=sin（x）/x
x2=[-0.57735,0.57735]
a2=[1,1];
x3=[-0.77460,0,0.77460]
a3=[0.55556,0.88889,0.55556]
x5=[-0.90618,-0.53847,0,0.53847,-0.90618]
a5=[0.23693,0.47863,0.56889,0.47863,0.23693]
fun1=[chebyshev（f1,2）,chebyshev（f1,3）,chebyshev（f1,5）]
fun2=[lanendre（f2,a2,x2）,lanendre（f2,a3,x3）,lanendre（f2,a5,x5）]
real1=double（int（f1/sqrt（1-x^2）,-1,1））
real2=double（int（f2,0,pi/2））
2点
3点
5点
真实值
Chebyshev
1.5708
Legendre
1.3704
1.3708
1.4327
1.3708
第九题
Euler法
function[t,x]=Euler（f,x0,h）
t=0:
h:
1;
m=length（t）;
x=zeros（1,m）;
x
（1）=x0;
m-1
x（i+1）=x（i）+h*subs（f,{'
t'
},{t（i）,x（i）}）;
改进Euler法
function[t,x]=ImproveEuler（f,x0,h）
k1=subs（f,{'
x1=x（i）+h*k1;
k2=subs（f,{'
},{t（i+1）,x1}）;
x（i+1）=x（i）+h/2*（k1+k2）;
Runge_Kutta法
function[t,x]=runge_kutta（f,x0,h）
k1=h*subs（f,{'
k2=h*subs（f,{'
},{t（i）+h/2,x（i）+k1/2}）;
k3=h*subs（f,{'
},{t（i）+h/2,x（i）+k2/2}）;
k4=h*subs（f,{'
},{t（i）+h,x（i）+k3}）;
x（i+1）=x（i）+（k1+2*k2+2*k3+k4）/6;
调用函数进行计算
symstx
f=x*cos（t）;
x0=1;
h=[0.1,0.01,0.001];
3
[t,x1]=Euler（f,x0,h（i））;
plot（t,x1）;
holdon;
>
clc
[t,x2]=Improved_Euler（f,x0,h（i））;
plot（t,x2）;
holdon;
[t,x3]=Runge_Kutta4（f,x0,h（i））;
plot（t,x3）;
此文档是由网络收集并进行重新排版整理.word可编辑版本！
升级会员 