数值计算方法 实验报告4

实验四数值微积分实验学院：数学与计算机科学学院专业：数学与应用数学学号：姓名：一．实验目的1利用复化求积公式计算定积分，并比较误差；2比较一阶导数和二阶导数的数值方法，并绘图观察特点.二．实验题目用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分，要求绝对误差为8105.0，并将计算结果与精度解进行比较：⑴dxexex2321432⑵dxxx322326ln.利用等距节点的函数值和端点的导数值，用不同的方法求下列函数的一阶和二阶导数，分析各种方法的有效性，并用绘图软件绘出函数的图形，观察其特点.⑴35611201xxy，2,0x⑵xey1，5.0,5.2x三．实验原理1复化梯形公式将积分区间ba,剖分为n等分，分点为)2,1,0(kkhaxk，其中nabh/)(.在每个区间1,kkxx上用梯形公式，则有dxxfdxxfnkxxbakk10110112nkkkkkkfRxfxfxxfRxfxfhnkknkkk101012.记10101222nkknkkknxfbfafhxfxfhT.2复化辛普森公式将积分区间ba,剖分为n等分，分点为)2,1,0(kkhaxk，其中nabh/)(.记区间1,kkxx的中点为21kx，在每个区间1,kkxx上用辛普森公式，则得到所谓的复化辛普森公式：12110146kkknkkknxfxfxfxxS，即102110426nkknkknxfxfbfafhS.3龙贝格公式的算法步骤为：I.输入ba,及精度；II.置,abhbfafhT211；III.置2,1,1nji，对分区间ba,，并计算111,ijijTT：nkkiixfhTT121111221，144111jijjjjijTTT；IV.若不满足终止条件，做循环：nnhhii2:,2/:,1:，计算nkkiixfhTT121111221，对,,,1ij计算：144111jijjjjijTTT.4向前差商公式：hafhafaf；向后差商公式：hhafafaf；中心差商公式：hhafhafaf2；二阶导数公式：22hhafafhafaf.四．实验内容实验一第一小题：对于方程dxexex2321432，利用程序shiyan1_01.m内容如下：%第一个函数的实验clearclcfun=inline('(2/3)*x.^3.*exp(x.^2)');S1=matrap(fun,1,2,170000);S2=masimp(fun,1,2,250);S3=maromb(fun,1,2,.5e-8);s=exp(4);Er1=abs(S1-s)Er2=abs(S2-s)Er3=abs(S3-s)第二小题：对于方程dxxx322326ln，利用程序shiyan1_02.m内容如下：%第二个函数的实验clearclcfun=inline('2*x./(x.^2-3)');S1=matrap(fun,2,3,15000);S2=masimp(fun,2,3,100);S3=maromb(fun,2,3,.5e-8);s=log(6);Er1=abs(S1-s)Er2=abs(S2-s)Er3=abs(S3-s)实验二第一小题：对于方程35611201xxy，2,0x，利用程序shiyan2_01.m内容如下：clearclcfun=inline('x.^5/20-(11./6)*x.^3');dfun=inline('x.^4/4-(11./2)*x.^2');ddfun=inline('x.^3-11*x');n=8;h=2/n;x=0:h:2;x1=x(2:n);y=feval(fun,x);dy=feval(dfun,x1);ddy=feval(ddfun,x1);fori=2:ndy1(i)=(y(i+1)-y(i))/h;dy2(i)=(y(i)-y(i-1))/h;dy3(i)=(y(i+1)-y(i-1))/(2*h);ddy1(i)=(y(i+1)-2*y(i)+y(i-1))/(h*h);endfori=1:n-1err1(i)=abs(dy1(i)-dy(i));err2(i)=abs(dy2(i)-dy(i));err3(i)=abs(dy3(i)-dy(i));errd2(i)=abs(ddy1(i)-ddy(i));end[err1'err2'err3'errd2']plot(x,y,'r')holdonplot(x1,dy,'y')plot(x1,ddy,'k')第二小题：对于方程xey1，5.0,5.2x，利用程序shiyan2_02.m内容如下：clearclcfun=inline('exp(-1./x)');dfun=inline('(-1./x).*exp(-1./x)');ddfun=inline('(-1./(x.^2)).*exp(-1./x)+1./(x.^2)');n=8;h=2/n;x=-2.5:h:-0.5;x1=x(2:n);y=feval(fun,x);dy=feval(dfun,x1);ddy=feval(ddfun,x1);fori=2:ndy1(i)=(y(i+1)-y(i))/h;dy2(i)=(y(i)-y(i-1))/h;dy3(i)=(y(i+1)-y(i-1))/(2*h);ddy1(i)=(y(i+1)-2*y(i)+y(i-1))/(h*h);endfori=1:n-1err1(i)=abs(dy1(i)-dy(i));err2(i)=abs(dy2(i)-dy(i));err3(i)=abs(dy3(i)-dy(i));errd2(i)=abs(ddy1(i)-ddy(i));end[err1'err2'err3'errd2']plot(x,y,'r')holdonplot(x1,dy,'y')plot(x1,ddy,'')五．实验结果实验一第一小题T=146.5012000000083.924363.065300000062.613255.509555.00580000056.653554.666954.610854.6045000055.115454.602754.598454.598254.598200054.727754.598454.598254.598254.598254.59820054.630554.598254.598254.598254.598254.598254.5982054.606254.598254.598254.598254.598254.598254.598254.5982Er1=4.5922e-009Er2=4.8409e-009Er3=1.4211e-014第二小题T=2.500000000002.01921.85900000001.85641.80221.7984000001.80881.79291.79221.792100001.79611.79181.79181.79181.79180001.79281.79181.79181.79181.79181.7918001.79201.79181.79181.79181.79181.79181.791801.79181.79181.79181.79181.79181.79181.79181.7918Er1=4.9383e-009Er2=4.0302e-009Er3=1.0132e-012实验二第一小题ans=0.21960.21960.21962.19200.36270.80030.58152.14800.57111.43671.00392.05600.76672.04111.40391.91600.94472.59911.77191.72801.10033.09632.09831.49201.22873.51832.37351.20801.32513.85072.58790.87601.38474.07912.73190.4960第二小题ans=0.69320.69320.69320.11050.46800.55320.51060.50300.52360.65550.58950.77930.59070.81020.70051.29910.66921.07270.87092.39820.74731.60711.17725.15720.75673.08731.922014.2888六．实验结果分析1.利用复化辛普森公式比利用复化梯形公式，所取的n更小，当达到相同精度时，利用辛普森公式等分次数n更小，减少计算次数.2.若利用同一公式，所取n的大小与题设给出的精度之间的关系：当n越大时，与精度之间的误差越小；反之，当n越小时，与精度之间的误差越大。3.龙贝格公式，是通过适当的线性组合，把复化梯形公式的近似值组合成更高精度的积分近似值.七．附件1复化梯形公式程序：%matrap.mfunctions=matrap(fun,a,b,n)%通途：用复化梯形公式求积分。%格式：s=matrap(fun,a,b,n)fun为被积函数，a,b为积分区间的左右端点，%n为区间的等分数，s返回数值积分值h=(b-a)/n;s=0;fork=1:n-1x=a+h*k;s=s+feval(fun,x);ends=h*(feval(fun,a)+feval(fun,b))/2.0+h*s;2复化辛普森公式程序：%masimp.mfunctions=masimp(fun,a,b,n)%通途：用复辛普生形公式求积分。%格式：s=masimp(fun,a,b,n)fun为被积函数，a,b为积分区间的左右端点，%n为区间的等分数，s返回数值积分值h=(b-a)/n;s1=0;s2=0;fork=1:(n-1)x=a+h*k;s1=s1+feval(fun,x);endfork=0:(n-1)x=a+h*(k+1/2);s2=s2+feval(fun,x);ends=h/6*(feval(fun,a)+feval(fun,b)+2*s1+4*s2);3龙贝格公式程序：%maromb.mfunctions=maromb(fun,a,b,tol)%用途：用龙贝格公式求积分%格式：s=maromb(fun,a,b,tol),fun是被积函数，a,b是积分下、上限%tol允许误差，s返回积分近似值ifnargin4,tol=1e-4;endi=1;j=1;h=b-a;T(1,1)=h*(feval(fun,a)+feval(fun,b))/2;T(i+1,1)=T(i,1)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2;T(i+1,j+1)=(4^j*T(i+1,j)-T(i,j))/(4^j-1);whileabs(T(i+1,i+1)-T(i,i))toli=i+1;h=h/2;T(i+1,1)=T(i,1)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2;forj=1:iT(i+1,j+1)=(4^j*T(i+1,j)-T(i,j))/(4^j-1);endendTs=T(i+1,j+1);