数值计算方法实验报告西安工程大学附有程序.docx
- 文档编号:3750605
- 上传时间:2022-11-25
- 格式:DOCX
- 页数:26
- 大小:100.16KB
数值计算方法实验报告西安工程大学附有程序.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)&(j j=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)’ 在根据Ux=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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 实验 报告 西安 工程 大学 附有 程序