数值计算方法实验报告西安工程大学附有程序.docx
《数值计算方法实验报告西安工程大学附有程序.docx》由会员分享，可在线阅读，更多相关《数值计算方法实验报告西安工程大学附有程序.docx（26页珍藏版）》请在冰豆网上搜索。
数值计算方法实验报告西安工程大学附有程序
数值分析
上
机
实
验
报
告
学院：
理学院
专业：
数学与应用数学
学号：
41108040110
姓名：
吕炳林
实验报告1
例1
实验程序
>>I0=exp（-1）*pond（‘x.^0.*exp（x.^2）’,0,1）
I0=
0.5381
>>vpa（I0,10）
实验结果：
ans=
0.5380795164
例1.9
实验程序：
functionp=p（n,x）
a
（1）=13;
fori=1:
n+1
a（i+1）=2*a（i）+3;
end
S=a（n+1）;
forj=1:
n
S=x*S+a（n+1-j）;
end
p=S;
End
在命令窗口进行计算：
实验结果
>>p（100,0.5）
ans=
600
>>p（150,13）
ans=
1.0995e+213
习题3
实验程序
Functiony=y（n）
Y=28
fori=1:
n
y=y-1/100*sqrt（783）;
End
在命令窗口进行计算
实验结果
y（100）
y=
28
ans=
0.0179
y（500）
y=
28
ans=
-111.9107
实验报告2
习题3
实验原理：
一次插值：
二次插值：
实验运行环境：
本实验采用Matlab编写。
实验程序：
线性插值：
x=0.4:
0.1:
0.8;
f=[-0.916291,-0.693147,-0.510826,-0.357765,-0.223144];
formatlong
interp1（x,f,0.54）
实验结果：
ans=
-0.620218600000000
二次插值：
x=0.54;
a=[0.4,0.5,0.6];
b=[-0.916291,-0.693147,-0.510826];
l=b
（1）*（x-a
（2））*（x-a（3））/（（a
（1）-a
（2））*（a
（1）-a（3）））;
m=b
（2）*（x-a
（1））*（x-a（3））/（（a
（2）-a
（1））*（a
（2）-a（3）））;
n=b（3）*（x-a
（1））*（x-a
（2））/（（a（3）-a
（1））*（a（3）-a
（2）））;
y=l+m+n
实验结果：
y=
-0.61531984000000
习题21
实验运行环境：
本实验采用Matlab编写
实验程序1：
（输入函数y）
forx=-4.5:
4.5
y=1/（x^2+1）
end
显示结果
y.0470********；
y=.0754********；
y=0.137********276；
y=0.30769230769231；
y=0.80000000000000；
实验程序2（求
的值）
x=input（'请输入x的值'）;
a=[x-0.5,x+0.5];
y=[1/（1+（x-0.5）^2）,1/（1+（x+0.5）^2）];
I=y
（1）*（x-a
（2））/（a
（1）-a
（2））+y
（2）*（x-a
（1））/（a
（2）-a
（1））
显示结果
当分别输入
时
I=0.0486I=0.0794I=0.1500I=0.3500I=0.7500
习题24
1
（2）
实验报告3
习题22
实验程序
x=[19,25,31,38,44]
y=[19.0,32.3,49.0,73.397.8]
[p,s]=polyfit（x,y,2）
实验结果p=
0.04970.01930.6882
s=
R:
[3x3double]
df:
2
normr:
0.1145
习题23。
实验程序
t=[0,0.9,1.9,3.0,3.9,5.0];
s=[0,10,30,50,80,110];
[P,S]=polyfit（t,s,5）
显示结果
P=
-0.54326.4647-26.560946.1436-13.2601-0.0000
S=
R:
[6x6double]
df:
0
normr:
1.2579e-012
实验结果：
习题24.
实验程序：
>>x=0:
5:
55;
y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64];
>>[P,S]=polyfit（x,y,5）
显示结果
P=
0.0000-0.00000.0002-0.00840.28510.0082
S=
R:
[6x6double]
df:
6
normr:
0.0487
实验结果：
实验报告4
习题4
实验程序
先求精确值：
x=sym（'x'）;
f=exp（-x）;
I=int（f,x,0,1）
I=
1-exp（-1）
利用simpson公式得；
a=0;b=1;
S=（b-a）/6*[exp（-a）+4*exp（-（a+b）/2）+exp（-b）]
显示结果S=0.6323
实验结果：
其误差为|S-I|=|0.6323+exp（-1）-1|≈1.7944e-004
习题8
实验程序；
ifnargin<5
Eps=1E-6;
end
m=1;
h=（b-a）;
err=1;
j=0;
T=zeros（4,4）;
T（1,1）=h*（limit（f,a）+limit（f,b））/2;
while（（err>Eps）&（jj=j+1;
h=h/2;
s=0;
forp=1:
m
x0=a+h*（2*p-1）;
s=s+limit（f,x0）;
end
T（j+1,1）=T（j,1）/2+h*s;
m=2*m;
fork=1:
j
T（j+1,k+1）=T（j+1,k）+（T（j+1,k）-T（j,k））/（4^k-1）;
end
err=abs（T（j,j）-T（j+1,k+1））;
end
I=T（j+1,j+1）;
ifnargout==1
T=[];
end
然后在输入框输入：
>>symsx;
>>f=sym（'exp（-x）*2/sqrt（3.1415926）'）
f=
exp（-x）*2/sqrt（3.1415926）
>>[I,T]=romberg（f,0,1,4,1E-8）
I=
0.7133
T=
0.77170000
0.72810.7135000
0.71700.71330.713300
0.71420.71330.71330.71330
0.71350.71330.71330.71330.7133
实验结果I=0.7133
习题11
Romberg方法
利用先前的Romberg.m文件
实验程序:
>>symsy;
f=sym（'1/y'）;
[I,T]=romberg（f,1,3,3,1E-5）
I=
1.0986
T=
1.33330000
1.16671.1111000
1.11671.10001.099300
1.10321.09871.09861.09860
1.09981.09861.09861.09861.0986
实验结果积分值为1.0986
（2）三点及五点Gauss公式
三点公式:
查表有,:
先作变换,将积分区间变换到上,令y=t+2,则有:
于是就可以利用求积公式来解题.
%三点公式
f=inline（'1/（t+2）'）
f=
Inlinefunction:
f（t）=1/（t+2）
>>x=[0.7745966692,-0.7745966692,0]
x=
0.7746-0.77460
>>A=[0.555555556,0.555555556,0.888888889]
A=
0.55560.55560.8889
>>I=0;
>>fori=1:
length（x）
I=I+feval（f,x（i））*A（i）;%三点Simpson公式
end
>>I
I=
1.0980
%五点公式
>>f=inline（'1/（t+2）'）
f=
Inlinefunction:
f（t）=1/（t+2）
>>x=[0.9061793,-0.9061793,0.5384693,-0.5384693,0]
x=
0.9062-0.90620.5385-0.53850
>>A=[0.2369263,0.2369263,0.4786287,0.4786287,0.5688889]
A=
0.23690.23690.47860.47860.5689
>>I=0;fori=1:
length（x）
I=I+feval（f,x（i））*A（i）;%三点Simpson公式
end
>>I
I=
1.0986
a=1;
b=3;
f=inline（'1/y'）;
%m=2,n=1Gauss-Legendre复化求积公式
n=1;
m=4;
h=（b-a）/（m）;
x=[0.5773502692-0.5773502692];
A=[11];
k=1;
I（k）=0;
forj=1:
n+1
temp=0;
fori=1:
m
tx=a+（2*（i-1）+1）*h/2+h*x（j）/2;
temp=temp+feval（f,tx）;
end
I（k）=I（k）+A（j）*temp;
end
I（k）=I（k）*h/2;
>>I
I=
1.0985
实验报告5
习题2实验程序：
（1）先用Euler方法：
>>F='x+y';
a=0;
b=1;
h=0.1;
n=（b-a）/h;
X=a:
h:
b;
Y=zeros（1,n+1）;
Y
（1）=0;
fori=2:
n+1
x=X（i-1）;
y=Y（i-1）;
Y（i）=Y（i-1）+eval（F）*h;
end
>>disp（[X',Y']）;
实验结果
ans=
01.0000
0.10001.1103
0.20001.2428
0.30001.3997
0.40001.5836
0.50001.7974
0.60002.0442
0.70002.3275
0.80002.6511
0.90003.0192
1.00003.4366
（2）用改进的Euler方法：
实验程序
>>F='x+y';
>>Y1=zeros（1,n+1）;
Y1
（1）=0;
fori=2:
n+1
x=X（i-1）;
y=Y1（i-1）;
ty=Y1（i-1）+eval（F）*h;
Y1（i）=Y1（i-1）+h/2*eval（F）;
x=X（i）;
y=ty;
Y1（i）=Y1（i）+h/2*eval（F）;
end
>>disp（[X',Y1']）;
实验结果ans=
01.0000
0.10001.1103
0.20001.2428
0.30001.3997
0.40001.5836
0.50001.7974
0.60002.0442
0.70002.3275
0.80002.6511
0.90003.0192
1.00003.4366
习题6
（1）实验程序
>>F='x+y';
a=0;
b=1;
h=0.2;
n=（b-a）/h;
X=a:
h:
b;
Y=zeros（1,n+1）;
Y
（1）=1;
fori=1:
n
x=X（i）;
y=Y（i）;
K1=h*eval（F）;
x=x+h/2;
y=y+K1/2;
K2=h*eval（F）;
x=x;
y=Y（i）+K2/2;
K3=h*eval（F）;
x=X（i）+h;
y=Y（i）+K3;
K4=h*eval（F）;
Y（i+1）=Y（i）+（K1+2*K2+2*K3+K4）/6;
end
实验结果；
>>disp（[X',Y']）;
01.0000
0.20001.2428
0.40001.5836
0.60002.0442
0.80002.6510
1.00003.4365
（2）
实验程序
>>F='3*y/（1+x）';
a=0;
b=1;
h=0.2;
n=（b-a）/h;
X=a:
h:
b;
Y=zeros（1,n+1）;
Y
（1）=1;
fori=1:
n
x=X（i）;
y=Y（i）;
K1=h*eval（F）;
x=x+h/2;
y=y+K1/2;
K2=h*eval（F）;
x=x;
y=Y（i）+K2/2;
K3=h*eval（F）;
x=X（i）+h;
y=Y（i）+K3;
K4=h*eval（F）;
Y（i+1）=Y（i）+（K1+2*K2+2*K3+K4）/6;
end
实验结果
>>disp（[X',Y']）;
01.0000
0.20001.7275
0.40002.7430
0.60004.0942
0.80005.8292
1.00007.9960
实验六
习题2
（1）实验程序
fa=’1-a*sin（a）’;
fc=’1-c*sin（c）’;
b=pi/2;
a=0;
R=1;
k=0；
while（R>5e-6）;
c=（a+b）/2;
ifeval（fa）*eval（fc）>0;
a=c;
else
b=c;
end
R=b-a;
k=k+1;
end
x=c
实验结果
ans=1.41021
习题4实验程序
fa=’exp（a）+10*a-2’;
fc=’exp（c）+10*c-2’;
b=1;
a=0;
R=1;
k=0；
while（R>5e-6）;
c=（a+b）/2;
ifeval（fa）*eval（fc）>0;
a=c;
else
b=c;
end
R=b-a;
k=k+1;
end
x=c
实验结果
x=0.090515136718750
（2）实验程序
f=input（'请输入需要求解函数>>','s'）
df=diff（f）;
miu=2;
x0=input（'inputinitialvaluex0>>'）;
k=0;max=100;R=eval（subs（f,'x0','x'））
while（abs（R）>1e-8）
x1=x0-miu*eval（subs（f,'x0','x'））/eval（subs（df,'x0','x'））;
R=x1-x0;
x0=x1;
k=k+1;
if（eval（subs（f,'x0','x'））<1e-10）;
break
end
ifk>max;
ss=input（'mayberesultiserror,chooseanewx0,y/n?
>>','s'）;
ifstrcmp（ss,'y'）
x0=input（'inputinitialvaluex0>>'）;
k=0;
else
break
end
end
end
k
x=x0
实验结果x=0.090526468052644
习题7实验程序
x=sym（'x'）;
f=sym（'x^3-3*x-1'）;
df=diff（f,x）;
FX=x-f/df
FX=x+（3*x-x^3+1）/（3*x^2-3）
Fx=inline（FX）;
x0=2;
fori=1:
10
disp（x0）;
x0=feval（Fx,x0）;
end
x0
formatshort;
显示结果
2
1.888888888888889
1.879451566951567
1.879385244836671
1.879385241571817
1.879385241571817
1.879385241571817
1.879385241571817
1.879385241571817
1.879385241571817
实验结果
x0=
1.879385241571817
习题13实验程序
f=input（'请输入需要求解函数>>','s'）
df=diff（f）;
miu=2;
x0=input（'inputinitialvaluex0>>'）;
k=0;max=100;R=eval（subs（f,'x0','x'））
while（abs（R）>1e-8）
x1=x0-miu*eval（subs（f,'x0','x'））/eval（subs（df,'x0','x'））;
R=x1-x0;
x0=x1;
k=k+1;
if（eval（subs（f,'x0','x'））<1e-10）;
break
end
ifk>max;
ss=input（'mayberesultiserror,chooseanewx0,y/n?
>>','s'）;
ifstrcmp（ss,'y'）
x0=input（'inputinitialvaluex0>>'）;
k=0;
else
break
end
end
end
k
x=x0
实验结果
所以
=1072380529
实验七
课题利用lu分解求出方程组CijXi=bi（其中i=j=100,C是一百阶方阵，b是从1到100的列向量）解线形方程组，根据题目要求C是一百阶方阵且满足下面约束条件
A=ones（100）;
C=diag（[diag（A,1）],-1）+diag（[diag（A,1）],1）+diag（[diag（A,2）],2）+diag（[diag（A,2）],-2）+4*eye（100）
b=（1:
100）’;
实验分析利用matlab中的lu函数对方阵C进行L,U分解
[L,U]=lu（C）
根据Ly=b可以求出y=inv（L）*b
实验结果
y=（1.0000,1.7500,2.4000,3.0000,3.6771,4.3475,5.0087,5.6727,6.3401,7.0060 , 7.6713 8.3372 ,9.0031 ,9.6689 ,10.3347, 11.0005 ,11.6663 ,12.3321 ,12.9980 13.6638 ,14.3296 ,14.9954 ,15.6612 ,16.3270 , 16.9928 ,17.6586, 18.3245 18.9903 ,19.6561 ,20.3219, 20.9877 ,21.6535 ,22.3193 ,22.9851, 23.6510 24.3168, 24.9826, 25.6484, 26.3142, 26.9800 ,27.6458, 28.3117 ,28.9775 29.6433, 30.3091, 30.9749 ,31.6407, 32.3065 ,32.9723 ,33.6382, 34.3040 34.9698, 35.6356 36.3014 36.9672 37.6330 38.2988 38.9647, 39.6305 40.2963, 40.9621, 41.6279, 42.2937, 42.9595 ,43.6254, 44.2912, 44.9570 45.6228, 46.2886 ,46.9544 ,47.6202 ,48.2860 , 48.9519, 49.6177, 50.2835 50.9493 ,51.6151 ,52.2809, 52.9467 ,53.6125, 54.2784, 54.9442, 55.6100 56.2758, 56.9416, 57.6074 ,58.2732 ,58.9390 ,59.6049, 60.2707 ,60.9365 61.6023, 62.2681, 62.9339 ,63.5997 ,64.2656, 64.9314, 65.5972, 66.2630 66.9288）’
在根据Ｕx=y可以求出x=inv（U）*y则方程的解为
实验结果
x=（0.0898, 0.2578, 0.3832 , 0.4960 , 0.6236 , 0.7514, 0.8751 ,0.9996， 1.1251 , 1.2501 , 1.3750 ,1.5000,1.6250 , 1.7500 , 1.8750 , 2.0000 2.1250 , 2.2500 ,2.3750 , 2.5000 , 2.6250 , 2.7500 , 2.8750 , 3.0000 3.1250 , 3.2500 , 3.3750, 3.5000 , 3.6250 , 3.7500 , 3.8750 , 4.0000 4.1250 , 4.2500, 4.3750 , 4.5000,4.6250 , 4.7500 , 4.8750 , 5.0000 5.1250 , 5.2500 , 5.3750 ,5.5000 , 5.6250 , 5.7500 ,5.8750 , 6.0000 6.1250 , 6.2500 , 6.3750 ,6.5000 , 6.6250 , 6.7500, 6.8750 , 7.0000 7.1250 ,7.2500 , 7.3750, 7.5000,7.6250 , 7.7500, 7.8750, 8.0000 8.1250 ,.2500, 8.3750 , 8.5000 , 8.6250 , 8.7500, 8.8750 , 9.0000, 9.1250 9.2500 ,9.3750 , 9.5000, 9.6250 ,9.7500, 9.8750, 10.0000, 10.1250 10.2501
升级会员 