计算方法实验报告

  • 格式:docx
  • 大小:68.64 KB
  • 文档页数:9

下载文档原格式

  / 9
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

计算方法实验报告(三)(一)数值微分法

一、实验问题

已知函数y=f(x)的数值如下,利用两点和三点数值微分公式计算和。

二、问题的分析(描述算法的步骤等)

基本的数值微分公式:

其中节点等距()且x为节点的特殊情形

两点数值微分公式:

=()/h

=()/h

三点数值微分公式:

=()/2h

二阶数值微分公式:

=()/

三、程序设计

两点数值微分公式——f0.m,f1.m

f0.m:

function y0 = f0( x1,x2,h )

y0=(x1-x2)/h;

end

f1.m:

function y1 = f1( x1,x2,h )

y1=(x1-x2)/h;

end

三点数值微分公式——f11.m

function y11 = f11( x2,x0,h )

y11=(x2-x0)/(2*h);

end

二阶数值微分公式——f2.m

function yy = f2( x2,x1,x0,h )

yy=(x2-2*x1+x0)/(h*h);

end

得出和结果的程序——weifen.m:

x=[0.6 0.8 0.9 1 1.1 1.2 1.4 ];

y=[0.7360 0.8365 0.9095 1 1.1105 1.2446 1.6017];

for i=1:3

h=0.1*2^(i-1);

d(1)=f0(y(4+i),y(4),h);

d(2)=f1(y(4),y(4-i),h);

d(3)=f11(y(4+i),y(4-i),h);

d(4)=f2(y(4+i),y(4),y(4-i),h);

fprintf('h=%.1f时:',h);

d=vpa(d,5)

end

四、计算结果

五、结果分析

随着h的增大,数值微分公式的截断误差越来越大。同时,在两点误差公式中很明显计算的误差与的误差的绝对值基本相等,同时在三点数值微分公式中求的误差在三个求一阶导数的公式中最小。

六、实验的总结与体会

数值微分公式的截断误差是随步长h的减小而减小的,但由于近似式本身均用h作除数,所以函数值计算的舍入误差,将随h的减小而给导数的计算带来越来越大的误差。因此h 的选取不能过小,原则上不能使舍入误差超过截断误差。在求解微分方程时,很明显三点数值微分方程结果更加精确。

(二)变步长积分法

一、实验问题

用变步长辛浦生积分公式和龙贝格积分法计算下列积分。

(1)(2)

二、问题的分析(描述算法的步骤等)

(1)变步长辛浦生积分法步骤,要求误差小于

L:令

对令

若则转,否则输出并停止计算

(2)龙贝格积分法步骤,要求误差小于:

对做

对令

对做

若D<则输出T并停止计算

打印失败信息,输出并停止计算三、程序设计

变步长辛浦生公式的函数——Simpson.m

function Svalue = Simpson( a,b,m )

h=(b-a)/2;

ss=f(a)+f(b);

sss=f(a+h);

ssss =0;

s=(ss+4*sss)*(h/3);

n=1;

while m>=n

h=h/2;

n=2*n;

ssss =sss+ssss;

如有帮助,欢迎下载支持。sss=0;

for i= 1:n;

sss=sss+f(a-h+2*i*h);

end

s=(ss+4*sss+2*ssss)*(h/3);

end

Svalue = vpa(s,7);

end

龙贝格积分法——Romberg.m

function Rvalue = Romberg (a,b,m)

h=b-a;

R(1,1)=(ff(a)+ff(b))*(h/2);

n=1;

j=0;

while m>=(n/2)

j=j+1;

h=h/2;

T=0;

for i=1:n

x=a+h*(2*i-1);

T=T+ff(x);

end

如有帮助,欢迎下载支持。

R(j+1,1)=R(j,1)/2+h*T;

for k=1:j

R(j+1,k+1)=R(j+1,k)+(R(j+1,k)-R(j,k))/(4^k-1);

end

n=2*n;

end

R=vpa(R,6)

Rvaule=vpa(R(j+1,j+1),6)

n

end

第(1)题中的函数——f.m

function y =f(x )

y=sin(x)/x;

if x == 0

y=1;

end

end

第(2)题中的函数——ff.m

function yy = ff( x )

yy=1/(1+x);

end

调用时的调用语句: