北科大工程数值上机实验.docx
《北科大工程数值上机实验.docx》由会员分享，可在线阅读，更多相关《北科大工程数值上机实验.docx（35页珍藏版）》请在冰点文库上搜索。
北科大工程数值上机实验
《工程数值计算》
上机实验报告
姓名丁恒
班级机自1404
学号41440001
2017年10月
1数据的插值与回归
1.1种插值算法的实践
1.待定系数法多项式插值
程序:
图像:
2拉格朗日多项式插值
x
0
1
2
y
1.0000
2.7183
7.3891
程序：
%拉格朗日线性插值函数的构造比较
clear;
x0=0;y0=1;
x1=1;y1=2.7183;
x2=2;y2=7.3891;
x=0:
0.01:
2;
fori=1:
length（x）
ye（i）=exp（x（i））;
L
（1）=（x1-x（i））*（x2-x（i））/（（x1-x0）*（x2-x0））;
L
（2）=（x0-x（i））*（x2-x（i））/（（x0-x1）*（x2-x1））;
L（3）=（x0-x（i））*（x1-x（i））/（（x0-x2）*（x1-x2））;
y（i）=y0*L
（1）+y1*L
（2）+y2*L（3）;
end
plot（x,y,'k-',x,ye,'r-',x0,y0,'ro',x1,y1,'ro'）;
结果：
3、牛顿多项式
程序：
clear;
xi=[0;40;80;120;160;200;240;280;320;360];
yi=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
n=length（xi）
B=zeros（n,n）;
B（1,:
）=yi';
fori=1:
n-1
forj=1:
n-i
B（i+1,j）=（B（i,j+1）-B（i,j））/（xi（i+j）-xi（j））;
end
end
x=1:
1:
360;
fori=1:
length（x）
ye（i）=log（x（i））;
y1（i）=B（1,1）+B（2,1）*（x（i）-xi
（1））;
y2（i）=B（1,1）+B（2,1）*（x（i）-xi
（1））+B（3,1）*（x（i）-xi
（1））*（x（i）-xi
（2））;
y3（i）=B（1,1）+B（2,1）*（x（i）-xi
（1））+B（3,1）*（x（i）-xi
（1））*（x（i）-xi
（2））+B（4,1）*（x（i）-xi
（1））*（x（i）-xi
（2））*（x（i）-xi（3））;
y4（i）=B（1,1）+B（2,1）*（x（i）-xi
（1））+B（3,1）*（x（i）-xi
（1））*（x（i）-xi
（2））+B（4,1）*（x（i）-xi
（1））*（x（i）-xi
（2））*（x（i）-xi（3））+B（5,1）*（x（i）-xi
（1））*（x（i）-xi
（2））*（x（i）-xi（3））*（x（i）-xi（4））;
end
plot（x,ye,'k-',x,y1,'r-',x,y2,'b-',x,y3,'g-',x,y4,'c-',xi,yi,'ro'）
结果：
n=10
4、分段线性插值
程序：
clear;
xi=[0;40;80;120;160;200;240;280;320;360];
yi=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
%主程序
n=length（xi）;
B=zeros（n,n）;
B（1,:
）=yi';
fori=1:
n-1
forj=1:
n-i
B（i+1,j）=（B（i,j+1）-B（i,j））/（xi（i+j）-xi（j））;
end
end
b=B（2:
end,1）;
%x=-5:
0.1:
5;
x=2:
0.1:
10;
forj=1:
length（x）
Y（1,j）=B（1,1）;
p（1,j）=1;
end
%牛顿插值
fori=1:
length（b）
forj=1:
length（x）
p（i+1,j）=p（i,j）*（x（j）-xi（i））;
Y（i+1,j）=Y（i,j）+b（i）*p（i+1,j）;
end
end
plot（xi,yi,'ko-',x,Y（4,:
）,'r-'）;
结果：
5、分段三次厄米特插值
插值程序：
clear;
N=10;
fori=1:
N+1
xi（i）=-5+10*（i-1）/N;
yi（i）=xi（i）^4;
ypi（i）=4*xi（i）^3;
end
n=length（xi）;
B=zeros（n,n）;
B（1,:
）=yi';
fori=1:
n-1
forj=1:
n-i
B（i+1,j）=（B（i,j+1）-B（i,j））/（xi（i+j）-xi（j））;
end
end
b=B（2:
end,1）;
x=-5:
0.1:
5;
forj=1:
length（x）
Y（1,j）=B（1,1）;
p（1,j）=1;
end
fori=1:
length（b）
forj=1:
length（x）
p（i+1,j）=p（i,j）*（x（j）-xi（i））;
Y（i+1,j）=Y（i,j）+b（i）*p（i+1,j）;
end
end
symsr
fori=1:
n-1
i
a=xi（i）;
b=xi（i+1）;
A（i,1）=（1+2*（r-a）/（b-a））*（（r-b）/（a-b））^2;
A（i,2）=（1+2*（r-b）/（a-b））*（（r-a）/（b-a））^2;
A（i,3）=（r-a）*（（r-b）/（a-b））^2;
A（i,4）=（r-b）*（（r-a）/（b-a））^2;
H（i,1）=A（i,1）*yi（i）+A（i,2）*yi（i+1）+A（i,3）*ypi（i）+A（i,4）*ypi（i+1）;
end
Hr=vpa（collect（H）,8）;
fori=1:
length（Hr）
Ce（i,:
）=sym2poly（Hr（i））;
end
fori=1:
length（x）
ye（i）=x（i）^4;
forj=1:
n-1
ifx（i）>=xi（j）&x（i）<=xi（j+1）
H3（i）=Ce（j,1）*x（i）^2+Ce（j,2）*x（i）+Ce（j,3）;
end
end
end
plot（x,ye,'k-',x,H3,'r--'）;
结果：
1/6分段三次样条插值
程序：
clear;
xi=[0;40;80;120;160;200;240;280;320;360];
yi=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
yp_L=0.5;
yp_R=-0.5;
n=length（xi）;
fori=1:
n-1
h（i）=xi（i+1）-xi（i）;
end
d=zeros（n,1）;
miu=zeros（n-2,1）;
lam=zeros（n-2,1）;
d
（1）=6/h
（1）*（（yi
（2）-yi
（1））/h
（1）-yp_L）;
d（n）=6/h（n-1）*（yp_R-（yi（n）-yi（n-1））/h（n-1））;
fori=2:
n-1
miu（i-1）=h（i-1）/（h（i-1）+h（i））;
lam（i-1）=1-miu（i-1）;
d（i）=6/（h（i-1）+h（i））*（（yi（i+1）-yi（i））/h（i）-（yi（i）-yi（i-1））/h（i-1））;
end
M=zeros（n,n）;
M（1,1）=2;
M（1,2）=1;
M（n,n-1）=1;
M（n,n）=2;
fori=2:
n-1
M（i,i）=2;
M（i,i-1）=miu（i-1）;
M（i,i+1）=lam（i-1）;
end
Sm=inv（M）*d;
symsr
fori=1:
n-1
S（i）=（xi（i+1）-r）^3/（6*h（i））*Sm（i）+（r-xi（i））^3/（6*h（i））*Sm（i+1）+（（yi（i+1）-yi（i））/h（i）-h（i）/6*（Sm（i+1）-Sm（i）））*（r-xi（i））+yi（i）-h（i）^2/6*Sm（i）;
end
Sr=vpa（collect（S）,8）;
fori=1:
length（Sr）
Se（i,:
）=sym2poly（Sr（i））;
end
x=0:
360;
fori=1:
length（x）
forj=1:
n-1
ifx（i）>=xi（j）&x（i）<=xi（j+1）
S3（i）=Se（j,1）*x（i）^3+Se（j,2）*x（i）^2+Se（j,3）*x（i）+Se（j,4）;
end
end
end
plot（xi,yi,'ko-',x,S3,'b-'）;
结果：
1.2MATLAB内部插值函数的应用
2.interp1函数
程序：
clear;
X=[0;40;80;120;160;200;240;280;320;360];
Y=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
xi=0:
360;
yi1=interp1（X,Y,xi,'*nearst'）;
yi2=interp1（X,Y,xi,'*linear'）;
yi3=interp1（X,Y,xi,'*spline'）;
yi4=interp1（X,Y,xi,'*cubic'）;
plot（X,Y,'ro',xi,yi1,'--',xi,yi2,'-',xi,yi3,'m:
'）
结果：
3.Spline函数
程序：
clear;
X=[0;40;80;120;160;200;240;280;320;360];
Y=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
xi=0:
360;
yi=spline（X,Y,xi）;
plot（X,Y,'o',xi,yi,'r-'）
结果：
4.Ppval函数
求20度、60度、100度处的函数值
程序：
clear;
xi=[0;40;80;120;160;200;240;280;320;360];
yi=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
p=csape（xi,yi,'variational'）;
x=[20,60,100];
ppval（p,x）
结果：
ans=
52.585144.269819.5631
1.3回归算法的实践
1、一维线性回归算法
问题：
x
1
2
3
4
y
2
2.5
3.33
4.25
程序：
%最小二乘法（线性拟合）
clear;
xi=[1;2;3;4];
yi=[2;2.5;3.33;4.25];
symsc1c2
n=length（xi）;
fori=1:
n
eb（i）=c1*xi（i）+c2-yi（i）;
end
J=0;
fori=1:
n
J=J+eb（i）^2;
end
J1=collect（J,c1）;
J2=collect（J,c2）;
diff（J,c1）
diff（J,c2）
x=1:
0.01:
4;
fori=1:
length（x）
y（i）=0.758*x（i）+1.125;
Y（i）=1/x（i）+x（i）;
end
plot（xi,yi,'ko',x,y,'r-',x,Y,'m-'）
结果:
1.4MATLAB内部回归函数的应用
1、polyfit和polyval函数
程序：
clear;
X=[0;40;80;120;160;200;240;280;320;360];
Y=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3];
p2=polyfit（X,Y,2）;
p3=polyfit（X,Y,3）;
p5=polyfit（X,Y,5）;
xi=[20,60,100,140];
yi2=polyval（p2,xi）
yi3=polyval（p3,xi）
yi5=polyval（p5,xi）
结果：
yi2=
47.610733.125722.745816.4708
yi3=
48.030831.806921.088515.4088
yi5=
56.554242.129819.82757.2451
2数值积分与数值微分
2.1数值积分算法的实践
1、梯形公式、辛普森公式和科特斯公式
程序：
clear;
xi=0:
5:
30;
F=[0.09.013.014.010.512.05.0];
du=[0.501.400.750.901.301.481.50];
f=F.*cos（du）;
C=zeros（3,4）;
C（1,1）=1/2;C（1,2）=1/2;
C（2,1）=1/6;C（2,2）=2/3;C（2,3）=1/6;
C（3,1）=1/8;C（3,2）=3/8;C（3,3）=3/8;C（3,4）=1/8;
a=0;
b=30;
%梯形公式
n=length（xi）;
T=0;
fori=1:
（n-1）
T=T+5*0.5*（f（i）+f（i+1））;
end
%辛普森公式
fori=1:
（n-2）
s（i）=C（2,1）*f（i）+C（2,2）*f（i+1）+C（2,3）*f（i+2）;
end
X=0;
fori=1:
（（n-1）/2）
X=X+10*s（2*i-1）;
end
%柯特斯公式
fori=1:
（n-3）
k（i）=C（3,1）*f（i）+C（3,2）*f（i+1）+C（3,3）*f（i+2）+C（3,4）*f（i+3）;
end
K=0;
fori=1:
（（n-1）/3）
K=K+15*k（3*i-2）;
end
结果：
T=
119.0892
>>X
X=
117.1271
>>K
K=
117.3265
2.2MATLAB内部数值积分函数的应用
分别用梯形公式和辛普森公式求函数y=4/（1+x^2）在[0,1]区间内的积分
1、quad函数、Trapz函数
y=4./（1+x.^2）
程序:
clear;
x=0:
0.01:
1;
y=4./（1+x.^2）;
%梯形公式
z1=trapz（x,y）;
%辛普森公式
z2=quad（'4./（1+x.^2）',0,1）;
z1=vpa（z1,8）
z2=vpa（z2,8）
结果：
z1=
3.141576
z2=
3.1415927
2.3数值微分算法的实践
1、对
求差商
程序：
不同步长对差商求导
clear;
h=zeros（17,1）;
h
（1）=0.05;
a0=1;
d0=1/（3*a0^（2/3））;
n=length（h）;
fori=2:
n
h（i）=h（i-1）/5;
end
fori=1:
n
a（i）=a0-h（i）;
b（i）=a0+h（i）;
c（i）=（b（i）^（1/3）-a（i）^（1/3））/（b（i）-a（i））;
e（i）=c（i）-d0;
er（i）=abs（e（i））/d0;
cp（i）=1/（a（i）^（2/3）+b（i）^（2/3）+（a（i）*b（i））^（1/3））;
ep（i）=cp（i）-d0;
erp（i）=abs（ep（i））/d0;
d（i）=d0;
end
subplot（2,1,1）;
semilogx（h,c,'r-',h,d,'ko-'）;
subplot（2,1,2）;
semilogx（h,er,'b-'）;
figure;
subplot（2,1,1）;
semilogx（h,cp,'r-',h,d,'ko-'）;
subplot（2,1,2）;
semilogx（h,erp,'b-'）;
结果：
2、对
求差商
程序：
%不同步长对差商求导
clear;
h=zeros（17,1）;
h
（1）=0.05;
a0=1;
d0=3*a0^2;
n=length（h）;
fori=2:
n
h（i）=h（i-1）/5;
end
fori=1:
n
a（i）=a0-h（i）;
b（i）=a0+h（i）;
c（i）=（b（i）^3-a（i）^3）/（b（i）-a（i））;
e（i）=c（i）-d0;
er（i）=abs（e（i））/d0;
cp（i）=a（i）^2+b（i）^2+a（i）*b（i）;
ep（i）=cp（i）-d0;
erp（i）=abs（ep（i））/d0;
d（i）=d0;
end
subplot（2,1,1）;
semilogx（h,c,'r-',h,d,'ko-'）;
subplot（2,1,2）;
semilogx（h,er,'b-'）;
figure;
subplot（2,1,1）;
semilogx（h,cp,'r-',h,d,'ko-'）;
subplot（2,1,2）;
semilogx（h,erp,'b-'）;
结果：
有上述结果知：
若步长h越小，截断误差越小，结果越精确
（1）当
，会由于截断误差而造成失真
（2）当
，则不会造成失真
2.4MATLAB内部数值微分函数的应用
diff函数的应用：
比较不同方法求得f（x）的导数
程序：
%数值微分函数的调用比较
x=-3:
0.1:
3;
f=inline（'sqrt（x.^3++2*x.^2-x+12）+（x+5）.^（1/6）+5*x+2'）;
g=inline（'（3*x.^2+4*x-1）./（2*sqrt（x.^3++2*x.^2-x+12））+1/6./（x+5）.^（5/6）+5'）;
p=polyfit（x,f（x）,5）;
dp=polyder（p）;
dpx=polyval（dp,x）;
dx=diff（f（[x,3.1]））/0.1;
gx=g（x）;
plot（x,g（x）,'r-',x,dpx,'bo',x,dx,'g--'）
结果：
3方程与方程组的求解
3.1线性方程组的求解
1、克莱姆法
程序：
%克莱姆法求解线性方程组
A=[223;477;-245];
b=[3;1;-7];
x=inv（A）*b
结果：
x=
2
-2
1
2、数值迭代法
程序：
%Jaboci迭代
A=[223;477;-245];
b=[3;1;-7];
n=length（b）;
N=200;
H=zeros（N,n）;
eb=1e-3
fori=1:
n
x（1,i）=0;
end
fori=1:
N
i
forj=1:
n
fork=1:
n
ifj~=k
H（i,j）=H（i,j）+A（j,k）*x（i,k）;
end
end
x（i+1,j）=（b（j）-H（i,j））/A（j,j）;
end
e（i）=abs（x（i+1,1）-x（i,1））;
ife（i）break
end
end
x
结果：
i=
20
x=000
1.50000.1429-1.4000
3.45710.6857-0.9143
2.1857-0.9184-0.5657
3.2669-0.54040.2090
1.7269-1.93290.3391
2.9243-1.18310.8371
1.4274-2.36530.7162
2.7910-1.38901.0632
1.2942-2.51520.8276
2.7738-1.42431.1298
1.2295-2.57200.8489
2.7986-1.40861.1494
1.1845-2.60580.8464
2.8362-1.38041.1584
1.1427-2.63630.8388
2.8781-1.34891.1661
1.0998-2.66790.8304
2.9223-1.31591.1742
1.0546-2.70120.8217
2.9687-1.28151.1829
结论：
数值迭代不收敛。
3.2非线性方程的求解
1、roots函数
程序：
%非线性方程的求根
p=[123];
roots（p）
结果：
ans=
-1.0000+1.4142i
-1.0000-1.4142i
2、fsolve函数
程序：
functionF=fun（x）
globalpr
F=p*sin（x）-r;
end
%非线性方程的求根
clear;
globalpr
p=3;
r=2;
x=fsolve（@fun,0）;
结果:
x=
0.7297
3、fzero函数
程序：
functiony=fun1（x）
y=exp（-x）-sin（pi/2*x）;
end%非线性方程的求根
clear;
x=fzero（@fun1,0）;
结果：
x=
0.4436
3.3非线性方程的求解
1、
程序：
%非线性方程组求解的应用
clear;
[x,y]=solve（'x^2+x*y+y-3=0','x^2-4*x+3=0'）
结果：
x=1
3
y=
1
-3/2
4微分方程的求解
初始条件：
程序：
sgm=10;
b=2.666667;
r=28;
x0=6;
y0=5;
z0=5;
T=20;
dt=0.005;
n=T/dt;
t=0:
dt:
20;
x
（1）=x0;
y
（1）=y0;
z
（1）=z0;
fori=1:
n
x（i+1）=x（i）+（-sgm*x（i）+sgm*y（i））*dt;
y（i+1）=y（i）+（r*x（i）-y（i）-x（i）*z（i））*dt;
z（i+1）=z（i）+（-b*z（i）+x（i）*y（i））*dt;
end
plot（x,y,'r-',x,z,'b-'）
结果：
、
初始条件：
结果:
4.1R-K算法及实践
程序：
%四阶R-K的应用
xi=0:
0.2:
1;
yi
（1）=1;
n=length（xi）;
h=（xi（n）-x
升级会员 