Romberg求积分公式

MATLAB程序设计实践课程评价1，为了实现以下科学计算算法的编程，同时以应用为例。Romberg积分公式2、编程解决了以下科学计算和工程实际问题：1)对于指定半径，为r，重量为q的统一圆柱体，轴的初始速度为v0，初始角速度为w0，v0r*w0，地面的摩擦系数为f，随着时间的推移，圆柱体无滑动滚动，从而得出此时圆柱轴的速度。2)，在丘陵上，高程在x和y方向每100米测量一次，高程数据如下：匹配曲面以确定合适的模型，并找到最高点和该点的高程。100200300400100636697624478200698712630478300680674598412400662626552334一、寻找Romberg积分公式1，算法说明：该算法预设了两个值的绝对或相对误差的容差。自动更改积分步长以使其小于Romberg加速方法公式对于等距离节点，可以通过自下而上除以球间距(a，b)的一半，然后乘以梯形公式来提高球坐标公式精度。这是一个幂级数形式，因为梯形公式的其馀部分只有精度，也就是说，节点加密时可以组合到达到其精度，与结合时可以达到误差精度，因此，取决于x，存在上阶导数的情况下可以扩展。利用外推原理，可以逐步消除右端的计算精度，通过将步骤h除以1/2，利用和组合消除，重复相同的过程，得到递归公式。此时，说明使用误差差来表示m次加速度。计算时以序列表示地块的一半数。也就是说，具体的计算公式是Romberg求积法。2，程序代码：m文件1)，Romberg加速法Function s，n=rbg1(a，b，eps)If nargin3，EPS=1e-6；EndS=10s0=0；k=2；T(1，1)=(B- a)*(f(a)f(b)/2；While (abs(s-s0)eps)h=(B- a)/2(k-1)；w=0；If (h=0)For I=1: (2 (k-1)-1)w=w f(a I * h)；EndT(k，1)=h *(f(a)/2 w f(b)/2)；For l=23360kFor i=1:(k-l 1)T (I，l)=(4 (l-1) * t (I 1，l-1)-t (I，l-1)/(4 (l-1)-)EndEndS=t(1，k)；S0=(t(1，k-1)；k=k 1；n=k；Else s=s0n=-k；EndEnd2)，改进的Romberg求积函数Function s，eer=rbg2(a，b，eps)If nargin3，EPS=1e-6；Endm=1；T(1，1)=(B- a)*(f(a)f(b)/2；R(1，1)=0；While (ABS (t (1，m)-r (1，m)/2) EPS)c=0；m=m 1；for j=1:2(m-1m-1)c=c f(a(j-0.5)*(B- a)/2(m-1)；EndR (m，1)=(B- a)* c/2(m-1)；For j=2:mFor k=1:(m-j 1)R (k，j)=r (k 1，j-1) (r (k 1，j-1)-r (k，j-1)/(4 (j-1)EndT (1，j)=r (1，j-1) 2 * (4 (j-2)-1) * (t (1，j-1)-r (1，)EndEndErr=abs(t(1，m)-r(1，m)/2；S=t(1，m)；3)将f.m函数定义为：function f=f(x)；f=x . 3；4)执行命令和结果Rbg1(0，2)Rbg2(0，2)3、流程图否是否是是否否是否是输入左端点和右端点a，b输入精度值EPSnargin 3Eps1e-6S=10s0=0；K=2T(1，1)=(b-a)*(f(a) f(b)/2Abs(s-s0)epsH=(b-a)/2 (k-1) w=0H=0？1=I=(2 (k-1)-1)W=w f(a i*h)T(k，1)=h*(f(a)/2 w f(b)/2)2=l=k？1=i=(k-l 1)否S=s0n=-k；流程图1I=i 1L=l 1I=12=L=2开始I=1T (I，l)=(4 (l-1) * t (I 1，l-1)-t (I，l-1)/(4 (l-1)-)S=t(k，1)；S0=(t(k-1，1)；k=k 1；n=k；导出结果结束I=i 1是是否是否是是否否是输入左端点和右端点a，b输入精度值EPSnargin 3Eps1e-6M=1T(1，1)=(b-a)*(f(a) f(b)/2(abs(t(1，m)-r(1，m)/2)epsc=0；m=m 1；J=11=j=2 (m-1)J=2否J=m流程图2J=j 1R(1，1)=0C=c f (a (j-0.5) * (b-a)/2 (m-1)J=j 1R (m，1)=(b-a) * c/2 (m-1)K=1K=m-j 1R (k，j)=r (k 1，j-1) (r (k 1，j-1)-r (k，j-1)/(4 (j-1)K=k 1T (1，j)=r (1，j-1) 2 * (4 (j-2)-1) * (t (1，j-1)-r (1，)导出结果结束Err=abs(t(1，m)-r(1，m)/2S=t(1，m)二、汽缸问题1、问题分析圆柱体的水平方向产生由地面的摩擦阻力f*Q加速到每个速度的力矩，该阻力充当轴速度的减速。总之，圆柱体没有摩擦，直到轴速度v=w*r。2，源：m文件r=input(r=)；Q=input(Q=)；g=input(g=)；f=input(f=)；v0=input(v0=)；w0=输入(w0=)；If v0 moveR=1Q=100G=9.81F=0.1V0=3W0=2T=0.3398V=2.6667三、高度4*4=有16个数据，每个数据对应16个坐标位置处的高程，您只需使用插值方法填充值即可拟合曲面。查找适合曲面的二进制函数更加复杂，并且插值后的数据量足够大(10000个)，因此，在求解山丘最高点及其高程时，可以比较所有数据，并将对应于该最大值的x，y值用作最高点。具体步骤如下2，源：m文件Function s，x0，y0=high(N)x=100 200 300 400；y=100 200 300 400；z=636 697 624 478；698 712 630 478；680 674 598 412；662 626 552 334；Xx=linspace(100，400，N)；Yy=linspace(100，400，N)；Zh=interp2(x、y、z、xx、yy和cubic)Mesh(xx，yy，zh)s=0；For i=1:N 2If zh(i)ss=zh(I)；N=mod(i，N)；m=(I-N)/N；EndEndX0=100 300/N*mY0=100 300/N*n3、执行命令高(100)运行结果如下：使拟合曲面如下求结果：Zh=Columns 1 through 7636 . 0000639 . 8173 643 . 5115 647 . 0826 650 . 5308 653 . 8558 657 . 0579639.0542 642.7968 646.4175 649.9164 653.2935 656.5488 659.6822642 . 0349 645 . 7037 649 . 2519 652 . 6794 655 . 9864659 . 1727 662 . 2384644 . 9421 648 . 5381652 . 0146 655 . 3717658 . 6094 661 . 727664 . 7265647 . 7759 654 . 7057 657 . 9932 661 . 1625 664 . 2136 667 . 1464650.5363 653.9892 657.3251 660.5440 663.6458 666.6305 669.4982653 . 2231656 . 6060 659 . 8729 663 . 0240 666 . 0591668 . 9784671 . 7818655 . 8365 659 . 1502 662 . 3491 665 . 4332668 . 4027 671 . 2574 673 . 9974658 . 3765 661 . 6218 664 . 7536 667 . 7717 670 . 6763673 . 4673 676 . 1448660 . 8430 664 . 0209 667 . 0864 670 . 0395672 . 8801 675 . 6083 678 . 2240663 . 2360 666 . 3475 669 . 3476 672 . 2365675 . 0140677 . 680 . 2352665 . 55668 . 6015 671 . 5372674 . 3627 677 . 0781679 . 6832 682 . 1781667 . 8017 670 . 7829 673 . 6551676 . 4182 679 . 0722 681 . 6172 684 . 0530669 . 9743672 . 8918675 . 7014 678 . 4029 680 . 9965 683 . 4821 685 . 8597672.0735 674.9282 677.6760 680.3169 682.8510 685.2781 687.5983674 . 0992 676 . 8920679 . 5790 682 . 1602 684 . 6355 687 . 0051 689 . 2688676.0514 678.7833 681.4103 683.9326 686.3502 688.6630 690.8711677.9302 680.6020 683.1700 685.6344 687.9950 690.2520 692.4053679 . 7355 682 . 3481 684 . 8581 687 . 2654 689 . 5700 691 . 7720 693 . 8713681 . 4674 684 . 0217 686 . 4745 688 . 8256 691 . 0751 693 . 2230695 . 2692683 . 1258 688 . 0192 690 . 3150 692 . 5103 694 . 6049 696 . 5990684 . 7107 687 . 1513 689 . 4923 691 . 7338 693 . 8756 695 . 9179 697 . 8607686 . 2222 688 . 6073 690 . 8938 693 . 0817695 . 1711 697 . 1619 699 . 0542687.6602 689.9907 692.2236 694.3589 696.3967 698.3369 700.1796689.0248 691.3016 693.4818 695.5654 697.5524 699.4429 701.2368690.3159 692.5399 694.6683 696.7011 698.6383 700.4799 702.2259691 . 5335 693 . 70576