数值分析实验报告1.docx
- 文档编号:6799959
- 上传时间:2023-01-10
- 格式:DOCX
- 页数:14
- 大小:54.50KB
数值分析实验报告1.docx
《数值分析实验报告1.docx》由会员分享,可在线阅读,更多相关《数值分析实验报告1.docx(14页珍藏版)》请在冰豆网上搜索。
数值分析实验报告1
数值分析上机实验报告
(注:
本实验报告中所有程序均为MATLAB语言程序)
班级:
姓名:
学号:
一章
1、利用数值积分计算
=
(n=0,1,2,……).
目的:
定积分数值求解
原理:
梯形公式法
程序:
clear
formatlong;
k=input('k=');
m=input('m=');
forn=1:
k
h=1/m;
x=0:
h:
1;
f=x.^n.*exp(x.^2);
fori=1:
m
s(i)=(f(i)+f(i+1))*h/2;
end
s=sum(s);
I(n)=exp(-1)*s;
end
I
运行结果:
k=9
m=1000
I=
Columns1through6
0.31606049880.23096057990.18394013730.15356013020.13212114220.1161015912
Columns7through9
0.10363907350.09364759740.0854476226
2、利用秦九韶算法计算当
=5,
=2
+3;n=100,x=0.5;n=150,x=13多项式
(
)=
+…
…
+
的值。
目的:
通过调整程序以简化计算步骤,减少运算次数
原理:
秦久韶算法
程序:
n=input('n=');
x=input('x=');
a
(1)=5;
fork=1:
n;
a(k+1)=2.*a(k)+3;
end
s(n+1)=a(n+1);
fori=n:
-1:
1
s(i)=x.*s(i+1)+a(i);
end
Pnx=s
(1)
运行结果:
n=100
x=0.5
Pnx=
802.0000000
n=150
x=13
Pnx=
1.4659714820e+213
3、设
=28,按递推公式
=
-1/100
计算到
,
。
若取
27.982(五位有效数字),试问计算
、
将有多大的误差。
目的:
使用递推公式计算求数列第n项值
原理:
用递推公式迭代
程序:
clear
k=input('k=')
Y0=28;
Y
(1)=Y0-sqrt(783)/100;
forn=2:
k
Y(n)=Y(n-1)-sqrt(783)/100;
end
y
(1)=Y0-27.982/100;
forn=2:
k
y(n)=y(n-1)-27.982/100;
end
y100=y(100)
y500=y(500)
e1=y(100)-Y(100)
e2=y(500)-Y(500)
运行结果:
k=500
k=
500
y100=
0.01800000000
y500=
-1.1191000000e+002
e1=
1.3715926637e-004
e2=
6.8579633158e-004
二章
目的:
掌握插值原理和不同插值方法
3、给出
=
的数值表,用线性插值级二次插值计算ln0.54的近似值。
x
0.4
0.5
0.6
0.7
0.8
Lnx
-0.916291
-0.693147
-0.510826
-0.357765
-0.223144
原理:
matlab线性插值函数
线性插值程序:
x1=input('x1=');
x=[0.4:
0.1:
0.8];
y=[-0.916291-0.693147-0.510826-0.357765-0.223144];
y1=interp1(x,y,x1,'linear')
运行结果:
x1=0.54
y1=
-0.620218600
原理:
二次插值
二次插值程序:
formatlong;
x0=[0.40.50.60.70.8];
y0=[-0.916291-0.693147-0.510826-0.356675-0.223144];
x=0.54;
n=length(x0);
s=0;
forj=0:
(n-1)
t=1;
fori=0:
(n-1)
ifi~=j
t=t*(x-x0(i+1))/(x0(j+1)-x0(i+1));
end
end
s=s+t*y0(j+1);
end
s
运行结果:
s=
-0.6161427152
4、给出cosx,
的函数表,步长h=
=
若函数表具有五位有效数字,研究用线性插值求cosx近似值时的总误差限。
原理:
线性插值
程序:
clear
i=1;
whilei<5402
x=0:
pi/180/60:
pi/2;
a=cos(x);
x1=x;
x(i)=[];
a(i)=[];
b(i)=interp1(x,a,x1(i),’linear’);
i=i+1;
end
b;
x=0:
pi/180/60:
pi/2;
a=cos(x);
a
(1)=[];
a(5400)=[];
b
(1)=[];
b(5400)=[];
c=b-a;
e1=max(abs(c))
运行结果:
e1=
4.2307972903e-008
21、求
=1/
在-5
5上取n=10,按等距节点求分段线性差值函数
,并计算各节点间终点处的
与
的值,并估计误差。
原理:
线性插值
程序:
formatlong
clear
x=-5:
5;
f=1./(1+x.^2);
fori=1:
10
Ihm(i)=(f(i)+f(i+1))/2;
xm(i)=(x(i)+x(i+1))/2;
fm(i)=1./(1+xm(i).^2);
e(i)=Ihm(i)-fm(i);
end
Ihm
fm
e
运行结果:
Ihm=
Columns1through6
0.048642533940.079411764710.15000000000.35000000000.75000000000.7500000000
Columns7through10
0.35000000000.15000000000.079411764710.04864253394
fm=
Columns1through6
.0470*******
Columns7through10
0.30769230770.13793103450.075471698110.04705882353
e=
Columns1through6
0.0015837104070.0039400665930.012068965520.04230769231-0.05000000000-0.05000000000
Columns7through10
0.042307692310.012068965520.0039400665930.001583710407
23.求f(x)=x^4在[a,b]上的分段Hermite插值,并估计误差。
原理:
matlab中Hermite插值函数
程序:
clear
a=input('a=');
b=input('b=');
n=input('n=');
x0=input('x0=');
x=linspace(a,b,n);
f=x.^4;
ifx0b
disp('Error');
else
y0=interp1(x,f,x0,'pchip');
end
y0
运行结果:
a=0
b=5
n=1000
x0=2
y0=
16.000000060
24,给定数据表如下,是求出三次样条插值S(x),并满足条件
原理:
三次样条差值
程序:
x=[0.250.300.390.450.53];
y=[0.50000.54770.62450.67080.7280];
x0=input('x0=');
y0=interp1(x,y,x0,'spline');(待定)
四章
4、用simpson公式求积分
并估计误差。
原理:
simpson公式法
程序:
clear
formatlong
a=0;
b=1;
f=inline('exp(-x)');
s=((b-a)/6)*(f(a)+4*f((a+b)/2)+f(b))
symsx
f=exp(-x);
s0=int(f,x,a,b);
s0=eval(s0);
e=s-s0
运行结果:
s=
0.6323336800
e=
2.131211751e-004
8、用Romberg方法计算积分2/
要求误差不超过
。
原理:
Romberg方法
首先建立一个以Romberg.m命名的M文件:
function[q,step]=Roberg(f,a,b,eps)
if(nargin==3)
eps=1.0e-4;
end;
M=1;
tol=10;
k=0;
T=zeros(1,1);
h=b-a;
T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b));
whiletol>eps
k=k+1;
h=h/2;
Q=0;
fori=1:
M
x=a+h*(2*i-1);
Q=Q+subs(sym(f),findsym(sym(f)),x);
end;
T(k+1,1)=T(k,1)/2+h*Q;
M=M*2;
forj=1:
k;
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
end;
tol=abs(T(k+1,j+1)-T(k,j));
end
q=T(k+1,k+1);
step=k;
运行结果:
>>f=@(x)exp(-x)
f=@(x)exp(-x)
>>[R,k,T]=romberg(f,0,1,1e-5)
R=0.6321
k=3
T=0.6839000
0.64520.632300
0.63540.63210.63210
0.63290.63210.63210.6321
Y=2./sqrt(pi).*R
Y=0.7133
11、用下列方法计算积分
,并比较结果。
(1)Romberg方法;
(2)三点及五点Gauss公式;
(3)将积分区间分为四等分,用复化两点Gauss公式。
Romberg方法:
借用上题建立的Romberg函数,在命令窗口输入:
>>f=@(y)1/(y)
f=@(y)1/(y)
>>[R,k,T]=romberg(f,1,3,1e-5)
R=1.098612291
k=5
三点Gauss公式
由于积分区间是[1,3],故作变化,
=
,查表有,
=0.77145966692,
。
程序:
f=inline('1/(x+2)');
x=[0.7745966692,-0.7745966692,0];
A=[0.555555556,0.555555556,0.888888889];
I1=0;
fori=1:
length(x)
I1=I1+feval(f,x(i))*A(i);
end
I1
运行结果:
I1=1.09803922
五点高斯公式:
=
,
=
,
,
=0.2369269,0.4786287,0.5688889
程序:
f=inline('1/(x+2)');
x=[-0.961799,-0.5384693,0,0.5384693,0.961799];
A=[0.2369269,0.4786287,0.5688889,0.4786287,0.2369269];
I2=0;
fori=1:
length(x)
I2=I2+feval(f,x(i))*A(i);
end
I2
运行结果:
I2=1.108682446
复化两点高斯公式:
同理
=
,通过查表求积系数为1,求积节点
。
程序:
f=inline('1/(x+2)');
I=1*feval(f,-0.57735)+1*feval(f,0.57735)
运行结果:
I=1.090908998
五章
未完待续……
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验 报告