利用Matlab实现Romberg数值积分算法----系统建模与仿真结课作业
使用Matlab进行数值积分的方法与注意事项
使用Matlab进行数值积分的方法与注意事项1. 引言数值积分是数学中的一个重要概念，它能够将曲线下的面积或者函数的总值进行估计和计算。
在实际应用中，由于很多函数无法直接进行解析求积，因此数值积分成为了一种常用的计算方法。
Matlab作为一款强大的数值计算软件，提供了很多用于数值积分的函数和方法。
2. 数值积分的基本原理数值积分的基本思想是将被积函数分割成一系列小区间，然后对每个小区间进行近似计算得到面积的总和。
这个过程可以看作是对大曲线的逼近，通过增加小区间的数目，可以得到更加精确的结果。
常见的数值积分方法有矩形法、梯形法、辛普森法等。
3. Matlab中的数值积分函数在Matlab中，有两个常用的数值积分函数分别是`quad`和`quadl`。
`quad`函数适用于一般的一元数值积分计算，而`quadl`函数则适用于具有奇点的积分计算。
这两个函数使用起来相对简单，只需要输入被积函数和积分区间即可。
例如，计算函数f(x)=x^2在区间[0, 1]上的积分可以使用以下代码：```f = @(x) x^2;integral = quad(f, 0, 1);disp(integral);```这段代码会输出函数f在区间[0, 1]上的积分值。
4. 数值积分的精度与误差控制在使用数值积分方法进行计算时，我们关心的一个重要问题是精度和误差控制。
数值积分的精度可以通过调整分割的区间数目来控制，一般来说，增加小区间的数目可以得到更加精确的结果。
此外，也可以通过提高数值积分方法的阶来提高精度。
Matlab中的`quad`和`quadl`函数具有较高的精度，并且可以通过设置选项来控制误差的允许范围。
5. 数值积分的注意事项在使用Matlab进行数值积分时，需要注意一些问题。
首先是积分区间的选择，需要确保被积函数在整个区间上是光滑的，没有奇点和间断。
如果存在奇点或者间断，需要通过分段积分或者奇点积分方法来处理。
其次是数值积分方法的选择，不同的函数可能适用于不同的数值积分方法，需要结合实际情况来选择最合适的方法。
用matlab实现romberg积分法
一、概述在数值分析中，求解定积分是一项重要的任务。
传统的数值积分方法包括梯形法则、辛普森法则等。
而Romberg积分法，是一种更加精确的数值积分方法，它通过不断增加区间的细分，逐步提高数值积分的精度。
在本文中，我们将尝试用MATLAB实现Romberg积分法，探索其优势和应用。
二、Romberg积分法原理Romberg积分法的基本原理是通过对梯形法则和辛普森法则进行逐步的细分和修正，以获得更加精确的数值积分结果。
假设我们需要求解函数 f(x) 在区间 [a, b] 上的定积分，那么Romberg积分法的步骤可以概括为以下几点：1. 将区间 [a, b] 均匀分成若干个子区间；2. 计算每个子区间上的梯形规则和辛普森规则的数值积分；3. 利用已知结果进行Richardson外推，修正数值积分的误差；4. 逐步增加子区间的细分，直到达到所需的精度要求。
三、MATLAB实现Romberg积分法我们可以使用MATLAB编程语言来实现Romberg积分法，以下是一个示例代码：function [I, R] = romberg(f, a, b, n)h = (b - a) ./ (2 .^ (0:n-1));R = zeros(n, n);R(1, 1) = (b - a) * (feval(f, a) + feval(f, b)) / 2;for j = 2:nsubtotal = 0;for i = 1:2^(j-2)subtotal = subtotal + feval(f, a + (2*i - 1) * h(j));endR(j, 1) = R(j-1, 1)/2 + h(j) * subtotal;endfor k = 2:nfor j = k:nR(j, k) = (4^(k-1) * R(j, k-1) - R(j-1, k-1)) / (4^(k-1) - 1); endendI = R(n, n);通过以上的MATLAB代码，我们可以轻松地实现Romberg积分法，对给定的函数和区间进行数值积分，并得到精确的积分结果。
数值积分算法与MATLAB实现
数值积分算法与MATLAB实现本文从网络收集而来，上传到平台为了帮到更多的人，如果您需要使用本文档，请点击下载按钮下载本文档（有偿下载），另外祝您生活愉快，工作顺利，万事如意！摘要：在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。
数值积分就是解决此类问题的一种行之有效的方法。
积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。
本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。
本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。
除了研究这些数值积分算法的理论外，本文还将这些数值积分算法在计算机上通过MATLAB软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。
【关键词】数值积分牛顿-科特斯求积公式高精度求积公式MATLAB软件前言对于定积分，在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式可以计算定积分的值，但在很多情况下的原函数不易求出或非常复杂。
被积函数的原函数很难用初等函数表达出来，例如等；有的函数的原函数存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式。
因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的。
另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值。
因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算。
而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值。
微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节。
数值分析实验— MATLAB实现
数值分析实验——MATLAB实现姓名sumnat学号2013326600000班级13级应用数学2班指导老师2016年1月一、插值：拉格朗日插值 (1)1、代码： (1)2、示例： (1)二、函数逼近：最佳平方逼近 (2)1、代码： (2)2、示例： (2)三、数值积分：非反常积分的Romberg算法 (3)1、代码： (3)2、示例： (4)四、数值微分：5点法 (5)1、代码： (5)2、示例： (6)五、常微分方程：四阶龙格库塔及Adams加速法 (6)1、代码：四阶龙格库塔 (6)2、示例： (7)3、代码：Adams加速法 (7)4、示例： (8)六、方程求根：Aitken 迭代 (8)1、代码： (8)2、示例： (9)七、线性方程组直接法：三角分解 (9)1、代码： (9)2、示例： (10)八、线性方程组迭代法：Jacobi法及G-S法 (11)1、代码：Jacobi法 (11)2、示例： (12)3、代码：G-S法 (12)4、示例： (12)九、矩阵的特征值及特征向量：幂法 (13)1、代码： (13)2、示例： (13)一、插值：拉格朗日插值1、代码：function z=LGIP(x,y)%拉格朗日插值n=size(x);n=n(2);%计算点的个数syms a;u=0;%拉格朗日多项式f=1;%插值基函数for i=1:nfor j=1:nif j==if=f;elsef=f*(a-x(j))/(x(i)-x(j));endendu=u+y(i)*f;f=1;endz=expand(u);%展开2、示例：>> x=1:6;y1=x.^5+3*x.^2-6;y2=sin(x)+sqrt(x);>> f1=LGIP(x,y1)f1 =-6+3*a^2+a^5%可知多项式吻合得很好>> f2=vpa(LGIP(x,y2),3)f2 =.962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5二、函数逼近：最佳平方逼近1、代码：function z=BestF(u,a,b,n)%最佳平方逼近，用x^i逼近,n为逼近的次数n=n+1;syms xreal;old=findsym(u);u=subs(u,old,x); %将u中变量替换为xf=sym('');H=sym('');d=sym('');for i=1:n %生成函数系f(1,i)=x^(i-1);endfor i=1:n %生成内积Hfor j=1:nH(i,j)=int(f(1,i)*f(1,j),a,b);endendfor i=1:n %生成内积dd(i,1)=int(f(1,i)*u,a,b);enda=H\d;%解法方程Ha=dz=a'*f';2、示例：>> syms x>> f1=sqrt(x);>> f2=x^5+x^2;>> f3=exp(x);>> a=0 ;b=1;>> BestF(f1,a,b,5)ans =12/143+420/143*x-1120/143*x^2+2016/143*x^3-1800/143*x^4+56/13*x^5>> BestF(f2,a,b,5)ans =x^5+x^2>> BestF(f3,a,b,5)ans =-566826+208524*exp(1)+(16733010-6155730*exp(1))*x+(-115830120+42611520*exp(1))* x^2+(306348840-112699440*exp(1))*x^3+(-342469260+125987400*exp(1))*x^4+(136302012-5 0142708*exp(1))*x^5>> vpa(ans,3)ans =.1e4-.1e6*x-.1e7*x^3+.1e7*x^4三、数值积分：非反常积分的Romberg算法1、代码：function z=IntRom(f,a,b) %Romberg 算法e=1e-10;I{1}=linspace(a,b,2);%1等分I{2}=linspace(a,b,3);%2等分L=setdiff(I{2},I{1});%新得插值点h=b-a;T(1,1)=h/2*sum(subs(f,I{1}));T(2,1)=1/2*T(1,1)+h/2*sum(subs(f,L));T(2,2)=4/3*T(2,1)-1/3*T(1,1);k=2;while abs(T(k,k)-T(k-1,k-1))>e %精度要求k=k+1;I{k}=linspace(a,b,2^(k-1)+1);L=setdiff(I{k},I{k-1});%集合差运算，新得插值点h=h/2;T(k,1)=1/2*T(k-1,1)+h/2*sum(subs(f,L));%梯形for i=2:kT(k,i)=(4^(i-1)/(4^(i-1)-1))*T(k,i-1)-(1/(4^(i-1)-1))*T(k-1,i-1);%加速endEndz=T(k,k);2、示例：>> syms x>> f=x^4;>> a=-4;b=4;>> IntRom(f,a,b)T =1.0e+003 *2.04800000000000 0 0 01.02400000000000 0.68266666666667 0 00.57600000000000 0.42666666666667 0.40960000000000 00.45200000000000 0.41066666666667 0.40960000000000 0.40960000000000ans =4.096000000000000e+002>> vpa((int(f,a,b)-ans),3)ans =0.>> f=exp(x);>> a=0;b=1;>> IntRom(f,a,b)T =Columns 1 through 41.85914091422952 0 0 01.75393109246483 1.71886115187659 0 01.72722190455752 1.71831884192175 1.71828268792476 01.72051859216430 1.71828415469990 1.71828184221844 1.718281828794531.71884112857999 1.71828197405189 1.71828182867536 1.718281828460391.71842166031633 1.71828183756177 1.71828182846243 1.71828182845905Columns 5 through 60 00 00 00 01.71828182845908 01.71828182845905 1.71828182845905ans =1.71828182845905>> vpa((int(f,a,b)-ans),3)ans =0.四、数值微分：5点法1、代码：function z=VDiff(f,x0)%5点法求导数值e=1e-15;h=0.01;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(1)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));%5点导数公式h=h/2;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(2)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));h=h/2;for i=-0:4x(i+1)=x0+i*h;endy=subs(f,x);m(3)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));k=3;while abs(m(k)-m(k-1))<abs(m(k-1)-m(k-2)) & abs(m(k)-m(k-1))>e & (h/10)>0%控制收敛条件及精度要求及h非0h=h/2;k=k+1;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(k)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));ende=abs(m(k)-m(k-1));z=[m(k);e];2、示例：>> syms x>> f=exp(x);>> x0=2;>> VDiff(f,x0)ans =7.389056098949710.00000000002558五、常微分方程：四阶龙格库塔及Adams加速法1、代码：四阶龙格库塔function z=RGFour(f,y0,a,b)%4阶龙格库塔，f为函数句柄h=0.01;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:n-1K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endz=Y;plot(X,Y);2、示例：>> f=@(x)sin(x(1));>> y0=0;a=0;b=2*pi;>> figure(1);>> RGFour(f,y0,a,b);3、代码：Adams加速法function z=myAdams(f,y0,a,b)h=0.01;p(4)=1;c(4)=1;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:3K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endfor i=4:n-1p(i+1)=Y(i)+h/24*(55*f([X(i),Y(i)])-59*f([X(i-1),Y(i-1)])+37*f([X(i-2 ),Y(i-2)])-9*f([X(i-3),Y(i-3)]));m(i+1)=p(i+1)+251/720*(c(i)-p(i));m(i+1)=f([X(i+1),m(i+1)]);c(i+1)=Y(i)+h/24*(9*f([X(i+1),m(i+1)])+19*f([X(i),Y(i)])-5*f([X(i-1), Y(i-1)])+f([X(i-2),Y(i-2)]));Y(i+1)=c(i+1)-19/720*(c(i+1)-p(i+1));endz=Y;plot(X,Y);4、示例：>> f=@(x)exp(x(1));>> myAdams(f,0,0,2*pi);六、方程求根：Aitken 迭代1、代码：function z=myAitken(f,x0);%Aitken 迭代求方程的根e=1e-15;xu1=x0;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);endwhile abs(xu2-xu1)>e%精度控制xu1=xu2;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);%Aitken加速公式endendz=xu2;2、示例：>> syms x>> f=cos(x/2)+x;>> x0=3;>> myAitken(f,x0)ans =3.14159265358979>> f=x^2-2+x;>> x0=1;>> myAitken(f,x0)ans =1.41421356237309七、线性方程组直接法：三角分解1、代码：function z=myGuess(A,b);%线性方程组三角分解求根n=size(A);if n~=rank(A)z=['矩阵A线性相关，不符合要求'];return;endn=n(2);L=eye(n);for i=1:n-1for j=i+1:nL(j,i)=A(j,i)/A(i,i);A(j,:)=A(j,:)-L(j,i)*A(i,:);endendU=A;for i=1:n %解Ly=b，得ys=0;for j=1:i-1s=s+y(j)*L(i,j);endy(i)=(b(i)-s)/L(i,i);endfor i=n:-1:1 %解Ux=y，得xs=0;for j=i+1:ns=s+x(j)*U(i,j);endx(i)=(y(i)-s)/U(i,i);endLUz=x';2、示例：>> A=[1 2 3;2 1 5;11 17 21];>> b=[1 3 5]';>> myGuess(A,b)L =1.00000000000000 0 02.00000000000000 1.00000000000000 011.00000000000000 1.66666666666667 1.00000000000000U =1.000000000000002.000000000000003.000000000000000 -3.00000000000000 -1.000000000000000 0 -10.33333333333333ans =-0.06451612903226-0.580645161290320.74193548387097>> t=A\bt =-0.06451612903226-0.580645161290320.74193548387097八、线性方程组迭代法：Jacobi法及G-S法1、代码：Jacobi法function z=myJacobi(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);%初始值e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=D\(L+U);f=D\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myJacobi(A,b)ans =3.000000000013402.000000000012610.999999999992373、代码：G-S法function z=myGS(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=(D-L)\U;f=(D-L)\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};4、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myGS(A,b)ans =3.000000000001351.999999999999160.99999999999954九、矩阵的特征值及特征向量：幂法1、代码：function z=myChar(A);%幂法求主特征值及对应的特征向量e=1e-10;n=size(A);n=n(2);v1=ones(n,1);u1=v1;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化while norm(u2-u1)>eu1=u2;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化endz{1}=c;z{2}=v2;2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> u=myChar(A)u =[14.00000000046956] [3x1 double]>> u{1}ans =14.00000000046956 >> u{2}ans =4.20000000191478 0.93333333198946 14.00000000046956。
数学实验“几种常见的求积分近似解的方法”实验报告（内含matlab程序）
数学实验“几种常见的求积分近似解的方法”实验报告（内含matlab程序）西京学院数学软件实验任务书实验二十一实验报告一、实验名称：Romberg 积分法，Gauss 型积分法，高斯-勒让德积分法，高斯-切比雪夫积分法，高斯-拉盖尔积分法，高斯-埃尔米特积分法。
二、实验目的：进一步熟悉Romberg 积分法，Gauss 型积分法，高斯-勒让德积分法，高斯-切比雪夫积分法，高斯-拉盖尔积分法，高斯-埃尔米特积分法。
三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言完成程序设计。
四、实验原理：1．Romberg 积分法：龙贝格积分法是用里查森外推算法来加快复合梯形求积公式的收敛速度，它的算法如下，其中()i m T 是通过一系列逼近原定积分的龙贝格分值.计算(0)1[()()]2b aT f a f b -=+ 对1,2,3,k n = ，计算下列各步：21()(1)11111(21)()[()]222k k k k k j b a j b a TT f a ---=---=++∑对1,2,,m k = 和,1,2,,1i k k k =-- ，计算111441m i i i m m m m T T T--+-=-随着计算的步骤的增加，()i mT 越来越逼近积分()ba f x dx ?。
2．Gauss 型积分法：高斯积分公式的思想是用n 个不等距的节点123,,,nx x x x 对被积函数进行插值，然后对插值后的函数进行积分，其积分公式为：111()()nk k k f x dx A f x -=≈∑?如果积分区间不是[1,1]-，则需转换到此区间：11()()222bab a b a b af x dx f t dt ---+=+?其中系数k A 、节点k x 与n 的关系如下表所示： 3．高斯-切比雪夫积分法：第一类切比雪夫积分形式为：11()()nk k k f x dx A f x -=≈∑?其中k A n π=，21cos2k k x nπ-= 4．高斯-拉盖尔积分法：高斯-拉盖尔公式有两种形式：1()()nxk k k e f x dx A f x +∞-=≈∑?1()()k nx k k k f x dx A e f x +∞=≈∑?下面编制的程序是针对第一种形式的高斯-拉盖尔公式，即1()()nxk k k e f x dx A f x +∞-=≈∑?因此程序的第一个输入参数——被积函数，是上式中的()f x 。
Romberg数值积分实验报告
Romberg数值积分实验报告班级：__ 姓名：学号：日期：_一、实验目的1、加深外推法的原理理解, 掌握Romberg外推法的计算方法。
2、用matlab软件实现Romberg数值积分来计算题目的运算。
二、基本理论及背景1、理论推导：对区间[a，b]，令h=b-a构造梯形值序列{T2K}。
T1=h[f(a)+f(b)]/2把区间二等分，每个小区间长度为h/2=(b-a)/2，于是T2 =T1/2+[h/2]f(a+h/2）把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是T4 =T2/2+[h/2][f（a+h/4）+f（a+3h/4).....................把[a，b] 2等分，分点xi=a+（b-a）/ 2 ·i （i =0，1，2 · · · 2k）每个小区间长度为（b-a）/ 2 .2、参考Romberg数值积分，实现积分的数值求解，完成下列题目：三、算法设计及实现1、算法设计(a)function y=fun1(x)y=sqrt(4-(sin(x)).^2)；(b)function y=fun2 (x)y=log(1+2*x)/3;(c)function y=fun3(x)y=exp(-x).*sin(x.^2)四、实验步骤1、○1打开matlab软件，新建Romberg.m文件,在窗口中编辑Romber数值积分函数程序代码,并保存在指定的文件夹下,在Current Directory窗口右边点击《Browse For Folder》按钮指向Romberg.m文件；○2在Command Window中编辑相应要计算的题目的数值函数及相应的题目的表达式。
2、输出结果和初步分析说明（见附件一）。
五、使用说明实验结果分析1、在Command Window窗口中编辑要调用的函数名与指定的函数名字不同导致出现错误，通过改正与函数名相同即可。
MATLAB数值计算绘图模拟仿真以及使用总结[全文5篇]
![MATLAB数值计算绘图模拟仿真以及使用总结[全文5篇]](https://img.taocdn.com/s3/m/8e2699e1370cba1aa8114431b90d6c85ec3a88dd.png)
MATLAB数值计算绘图模拟仿真以及使用总结[全文5篇]第一篇：MATLAB数值计算绘图模拟仿真以及使用总结Work1 1-11-21-31-4Work2 2-12-1-（1）2-22-32-3-(1)2-4 2-4-(1)2-52-6和2-7Work3 3-13-1-（1）Work 4 4-1 4-24-2-（2）4-2-3Work5 5-1-（1）5-1-25-1-（3）5-2简述MATLAB命令窗的主要作用？（1）命令窗口(Command Window)位于MATLAB 操作桌面的右方，用于输入命令并显示除图形以外的所有执行结果，是MATLAB 的主要交互窗口。
（2）Matlab既可以运行命令也可以执行程序，在命令窗口中可以运行单独的命令也可以调用程序，相当方便，而编辑调试窗口和图像窗口都是程序运行结果展示窗口，可以很直观的对程序运行过程中出现的矩阵或者是变量等等进行监视。
（3）在MATLAB 命令窗口中可以看到有一个“>>”，该符号为命令提示符，表示MATLAB正在处于准备状态。
在命令提示符后输入命令并按回车键后，MATLAB 就会解释执行所输入的命令，并在命令后面给出计算结果。
5-3简述MATLAB绘制二维图形的一般步骤 MATLAB绘制图形一般采取以下7个步骤：（1）准备数据（2）设置当前绘图区（3）绘制图形（4）设置图形中曲线和标记点格式（5）设置坐标轴和网格线（6）标注图形（7）保存和导出图形5-4启动Simulink的方式有几种？ 1.启动Simulink 启动Simulink通常有三种方式：1)直接从Matlab指令窗口选取菜单File| New| Modal命令，Matlab将会打开Simulink库浏览器和名为untitled的模型窗口。
2)在Matlab命令窗口中键人Simulink命令，Matlab将会打开Simulink 库浏览器。
3)点击Matlab命令窗口工具条的图标，启动Simulink库浏览器。
基于matlab的蒙特卡罗积分的实现
概率论与数理统计论文基于matlab的蒙特卡罗定积分的实现班级：姓名：基于matlab的蒙特卡罗定积分的实现摘要：在对蒙特卡罗方法概念学习的基础上，利用计算机产生了随机数序列．在此基础上，研究了蒙特卡罗积分，并应用matlab加以实现．关键词：蒙特卡罗；MATLAB；定积分0 引言随着电子计算的出现和发展，近年来用概率模型来作近似计算的方法得到了很大的发展，即蒙特卡罗(Monte—Garlo)方法．它是一种采用统计抽样理论近似的求解数学与物理问题的方法，它既可以用来研究概率问题，也可以用来解决非概率问题．蒙特卡罗法已被广泛地运用到各个领域中，如高维数学问题求解、医学技术中的诊断识别、大型系统的可靠性分析等．一般蒙特卡罗方法在数学中最常见的应用就是蒙特·卡罗积分．对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法能够提供一种有效而可行的解决方法．1 蒙特卡罗方法的基本原理用蒙特卡罗方法解决数学分析问题时，基本思想是：首先建立与描述该问题有相似性的概率模型，利用这些相似把使该模型的一些数字特征(如概率或均值)与数学分析的解答(如积分值或微分方程的解)联系起来；然后对该模型进行随机模拟或统计抽样；最后利用所得的结果求出这些特征的统计估计值作为原来问题的近似解．蒙特卡罗方法计算的结果收敛的理论依据来自于大数定律，且结果渐进地服从正态分布的理论依据是中心极限定理．本文将对蒙特卡罗方法计算定积分，并用matlab 加以实现．2 蒙特卡罗方法求定积分的一个实例2．1 采用均匀随机数的蒙特卡罗方法计算和 matlab 实现 为了使用蒙特卡罗方法求该积分dx x I ⎰=13，首先选定一个矩形区域 D =[0，1；O ，1]．然后抽取n 个随即点的坐标()i i ηξ,(i=1，2，⋯，n)，他们服从区域D 上的二维均匀分布．计算落在区域D 内的随机点数设为m ，即满足0≤i η≤3x．则由强大数定律，要求的定积分近似的等于两个区域随机落点的比值 与矩形面积的乘积． matlab 程序及其运行结果： 2．2 matlab 生成随机掷点效果图 x1=rand(1,2000); y1=rand(1,length(x1)); x=0:0.01:1; y2=0;y3=x.^3;subplot(2,1,1),plot(x,y2,'k',x,y3,'k',x1,y1,'.') xlabel('x 轴'),ylabel('y 轴') m=1;for j=10:10:30000;a(m)=j;x0=rand(1,j);y0=rand(1,length(x0));s=0;n=length(x0);for i=1:n y4(i)=0;y5(i)=x0(i).^3; if y0(i)<y5(i)&y0(i)>y4(i) s=s+1; endendp(m)=s/n;m=m+1;endt=mean(p);txx=0;length(x0);yy=1/4;subplot(2,1,2),plot(a,p,'b',xx,yy,'b')t ≈ 0.24992．3 matlab 均匀随机数法计算定积分由以上运行结果可知，均匀随机数法求的定积分dx x I ⎰=13≈ 0.2499由定积分的几何意义可知，它是三次曲线3x y =与直线y = 0所围的平面图形的面积，由微积分基本公式可知，4141410103===⎰x dx x I ．模拟值与精确值的误差率为0.4％．3 小结Monte Carlo 方法的实质是通过大量随机试验，利用概率论解决问题的一种数值方法，基本思想是基于概率和体积间的相似性，因此程序十分简单；蒙特卡罗方法在解决确定性问题时，具有较高的精确度；误差主要取决于样本的容量而与维数无关，随着抽取点增加，近似面积也将逼近真实面积；对问题求解的过程取决于所建立的概率模型．参 考 文 献：[1] 概率论与数理统计教程/茆诗松，程依明，濮晓龙编著. —2版. —北京：高等教育出版社，2011.2[2] MATLAB 教程及实训/曹弋主编. —北京：机械工业出版社，2008.4。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
利用Matlab 实现Romberg 数值积分算法
一、内容摘要
针对于某些多项式积分，利用Newton —Leibniz 积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab 中的.m 文件编写了复化梯形公式与Romberg 的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。
二、数值积分公式
1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton —Cotes
求积公式:
I =(x)[f(a)f(b)]2
b
a b a
f dx -≈
+⎰ 其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f 在区间[a,b]上的积分值，截断误差为：
3"
(b a)()12
f η-- (a,b)η∈ 具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：
I =1
1
(b a)(b a)
(x)dx [f(a)f(b)2(a )]2n b
a
k k f f n n -=--≈+++∑⎰
其截断误差为：
2"
(b a)h ()12
R f η--=
(a,b)η∈ 2.Romberg 数值积分算法
使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg 数值积分。
其思想主要是，根据I 的近似值2n T 加上I 与2n T 的近似误差，作为新的I 的近视，反复迭代，求出满足计算精度的近似解。
用2n T 近似I 所产生的误差可用下式进行估算：
12221
()3
n n n I T T T -∆=-=-
新的I 的近似值：
122
n n j T T -=∆+ j =（0 1 2 ….） Romberg 数值积分算法计算顺序
i=0 (1) 002T
i=1 (2) 102T
(3) 012T
i=2 (4) 202T (5) 112T
(6) 022T
i=3
(7)
302T (8) 212T (9) 122T
(10) 032T i=4 (11)
402T
(12) 312T
(13) 222T
(14) 132T
…
…
…
…
其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg 序列。
三、复化梯形法以及Romberg算法程序流程图
图1 复化梯形法程序流程图
图2 Romberg算法程序流程图
四、计算实例
依据上文所述的流程图，编写复化梯形程序以及Romberg 算法程序，并且利用实例验证程序的正确性，示例如下(计算精度)：
1
2
04
1dx x π=+⎰
表2 计算结果
计算精度 0.5×10^-5 0.5×10^-7
0.5×10^-9
复化梯形 算法 时间
0.069826394633 0.216635802304 3.459824945493
近似值 3.141590110458 3.141592613853 3.141592653434 Romberg 算法
时间
0.0456******** 0.0433******** 0.044913907518
近似值 3.141592502458 3.141592651224 3.141592653552
从上表中可以看出，当要求的计算精度不高时，复化梯形算法与Romberg 算法计算时间相差不太大，但是Romberg 算法是要快于复化梯形算法的；当要求的计算精度更高的时候，Romberg 算法是明显快于复化梯形算法。
本文所编写的程序适用于多项式的数值积分，且对于积分区间内，被积函数在每一点必须有定义，在以后的学习中进一步改进。
附录：
1.复化梯形算法程序
function []=sf(a,b,m,M,d)
tic
disp('请输入分子多项式a，分母多项式b，积分下限m,积分上限M,以及计算精度d')
f=poly2sym(a)/poly2sym(b) %用于给用户显示被积函数的形式%利用梯形公式计算此数值积分
disp('利用梯形公式计算数值积分的结果')
kk=zeros(); %用于存放结果
kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M)) %先存储首项
for i=1:1:2^30
t=0;
for j=0:1:2^(i-1)-1
v=m+(2*j+1)*(M-m)/(2^i)
vv=polyval(a,v)/polyval(b,v);
t=t+(M-m)/(2^i)*vv
end
y=1/2*kk(i,1)+t %通项公式计算各项值
kk(i+1,1)=y %存储其他项
f=i+1; %记录符合条件的值的下标
if(1/3*(kk(i+1,1)-kk(i,1))<=d)
break;
end
end
time=toc
fprintf('The result is %f\n', kk(f,1))
2.Romberg算法程序
function []=romberg(a,b,m,M,d)
tic
disp('请输入分子多项式a，分母多项式b，积分下限m,积分上限M,以及计算精度d')
f=poly2sym(a)/poly2sym(b) %用于给用户显示被积函数的形式disp('利用梯形公式计算数值积分的结果')
kk=zeros(); %用于存放结果
kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M)); %先存储首项
for i=1:1:2^40
t=0;
for j=0:1:2^(i-1)-1
v=m+(2*j+1)*(M-m)/(2^i);
vv=polyval(a,v)/polyval(b,v);
t=t+(M-m)/(2^i)*vv;
end
y=1/2*kk(i,1)+t; %通项公式计算各项值
kk(i+1,1)=y ; %存储其他项
if(abs(1/3*(kk(i+1,1)-kk(i,1)))<=d) %判断梯形公式值是否达到要求disp('The result is:')
kk()
kk(i+1,1) %梯形值满足要求，输出结果
break;
else
s=(4*kk(i+1,1)-kk(i,1))/(4-1); %构造simpson各项
kk(i+1,2)=s %存储
if(i+1>=3)
if(i+1>=3 & abs(1/15*(kk(i+1,2)-kk(i,2)))<=d)
kk()
disp('The result is:')
kk(i+1,2) %simpson值满足要求，输出结果
pan1=0;
break;
else
c=(4^2*kk(i+1,2)-kk(i,2))/(4^2-1);%构造cotes值
kk(i+1,3)=c %存储cotes值
if(i+1>=4)
if(i+1>=4 & abs(1/63*(kk(i+1,3)-kk(i,3)))<=d)
disp('The result is:')
kk(i+1,3)
break;
else
r=(4^3*kk(i+1,3)-kk(i,3))/(4^3-1)%构造romberg值
kk(i+1,4)=r %存储romberg值
if(i+1>=5)
if(i+1>=5 & abs(1/127*( kk(i+1,4)- kk(i,4)))<=d )
disp('The result is:')
kk(i+1,4)
break;
end
end
end
end
end
end
end
end
time=toc。