利用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。
如何使用Matlab进行系统建模和仿真
如何使用Matlab进行系统建模和仿真一、引言在现代科学和工程领域，系统建模和仿真是解决实际问题和优化设计的重要手段之一。
Matlab作为一种功能强大的工具，被广泛应用于系统建模和仿真。
本文将介绍如何使用Matlab进行系统建模和仿真的基本步骤，并通过实例演示其应用。
二、系统建模系统建模是将实际系统抽象成数学或逻辑模型的过程。
在Matlab中，可以使用符号表达式或差分方程等方式对系统进行建模。
1. 符号表达式建模符号表达式建模是一种基于符号计算的方法，可以方便地处理复杂的数学运算。
在Matlab中，可以使用符号工具箱来进行符号表达式建模。
以下是一个简单的例子：```matlabsyms xy = 2*x + 1;```在上述例子中，定义了一个符号变量x，并使用符号表达式2*x + 1建立了y的表达式。
通过符号工具箱提供的函数，可以对y进行求导、积分等操作，从而分析系统的特性。
2. 差分方程建模差分方程建模是一种基于离散时间的建模方法，适用于描述离散时间系统。
在Matlab中，可以使用差分方程来描述系统的行为。
以下是一个简单的例子：```matlabn = 0:10;x = sin(n);y = filter([1 -0.5], 1, x);```在上述例子中，定义了一个离散时间信号x，通过filter函数可以求得系统响应y，其中[1 -0.5]表示系统的差分方程系数。
三、系统仿真系统仿真是利用计算机模拟系统的运行过程，通过数值计算得到系统的输出响应。
在Matlab中，可以使用Simulink工具箱进行系统仿真。
1. 搭建系统框图在Simulink中，我们可以使用各种模块来搭建系统的框图。
例如，可以使用连续时间积分器模块和乘法器模块来构建一个简单的比例积分控制器：在上图中，积分器模块表示对输入信号积分，乘法器模块表示对输入信号进行放大。
利用Matlab实现Romberg数值积分算法----系统建模与仿真结课作业
利用Matlab 实现Romberg 数值积分算法一、内容摘要针对于某些多项式积分，利用Newton —Leibniz 积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab 中的.m 文件编写了复化梯形公式与Romberg 的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。
二、数值积分公式1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton —Cotes求积公式:I =(x)[f(a)f(b)]2bab af dx -≈+⎰ 其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f 在区间[a,b]上的积分值，截断误差为：3"(b a)()12f η-- (a,b)η∈ 具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：I =11(b a)(b a)(x)dx [f(a)f(b)2(a )]2n bak k f f n n -=--≈+++∑⎰其截断误差为：2"(b a)h ()12R f η--=(a,b)η∈ 数值积分算法使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg 数值积分。
其思想主要是，根据I 的近似值2n T 加上I 与2n T 的近似误差，作为新的I 的近视，反复迭代，求出满足计算精度的近似解。
用2n T 近似I 所产生的误差可用下式进行估算：12221()3n n n I T T T -∆=-=-新的I 的近似值：122n n j T T -=∆+ j =（0 1 2 ….） Romberg 数值积分算法计算顺序i=0 (1) 002Ti=1 (2) 102T(3) 012Ti=2 (4) 202T (5) 112T(6) 022Ti=3(7)302T (8) 212T (9) 122T(10) 032T i=4 (11)402T(12) 312T(13) 222T(14) 132T…………其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg 序列。
数值积分的matlab实现
然后编写如下程序,并保存为t2.m
%t2.m
clear;
a=0;b=2*pi;
T=chen_trap(@f2,a,b);
S=chen_simpson(@f2,a,b);
FT8=comtrap(@f2,a,b,8);
FS4=comsimpson(@f2,a,b,8);
quad=Romberg(@f2,a,b,10e-6);
(4)
(5)
1、依照实验原理编写各种公式的的程序.
2、首先在电脑上安装matlab，然后启动matlab，分别建立不同公式的M文件，
实验程序如下；
梯形公式程序代码
程序代码说明
function T=chen_trap
(fun,a,b)；
T=(b-a)*(feval(fun,a)+feval(fun,b))/2；
disp('与精确值比较的误差为:'),
wuchi=[double(abs(Sj-T)),double(abs(Sj-S)),double(abs(Sj-FT8)),double(abs(Sj-FS4)),double(abs(Sj-quad))]
五、实验数据及结果：
运行t1得
t1
积分准确值为:
Sj = 0.492712
disp('与精确值比较的误差为:'),
wuchi=[double(abs(Sj-T)),double(abs(Sj-S)),double(abs(Sj-FT8)),double(abs(Sj-FS4)),double(abs(Sj-quad))]
;
首先编写被积函数的M文件f2.m
function y=f2(x)
(2021年整理)数值积分用matlab实现
数值积分用matlab实现编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（数值积分用matlab实现）的内容能够给您的工作和学习带来便利。
同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。
本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快业绩进步，以下为数值积分用matlab实现的全部内容。
东北大学秦皇岛分校数值计算课程设计报告数值积分及Matlab实现学院数学与统计学院专业信息与计算科学学号5133117姓名楚文玉指导教师张建波姜玉山成绩教师评语：指导教师签字：2015年07月14日1 绪论在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题． 1.1 课题的背景对于定积分()ba f x dx ⎰在求某函数的定积分时,在一定条件下，虽然有牛顿—莱布里茨公式()()()ba I f x dx Fb F a ==-⎰可以计算定积分的值，但在很多情况下的原函数()f x 不易求出或非常复杂．被积函数的原函数很难用初等函数表达出来,例如2sin (),x x f x e x-=等;有的函数()f x 的原函数()F x 存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式．因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的．另外,许多实际问题中的被积函数()f x 往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值．因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算．而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值．微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一．随着计算机的出现，近几十年来,对于数值积分问题的研究已经成为一个很活跃的研究领域．现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义．国内外众多学者在数值积分应用领域也提出了许多新方法．在很多实际应用中,只能知道积分函数在某些特定点的取值,比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等．通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab 软件编程的实现，应用于实际生活中． 1.2 课题的主要内容框架1。
龙贝格数值积分报告(MATLAB)
目录设计总说明 (1)关键词 (1)绪论 (2)第1章课程设计正文 (3)1.1 需求分析 (3)1.2 界面设计 (3)1.3 功能实现 (4)1.4系统算法和程序代码的实现 (5)第2章总结 (7)参考文献 (8)设计总说明学习了《数值分析》和相关开发软件课程并能够熟练掌握各种算法的原理及其求解方法和特点；1、独立进行面向对象软件的分析、设计及在计算机中进行编写代码和调试工作。
培养自己独立进行面向对象软件的分析、面向对象软件的设计及面向对象软件的实现与测试的能力；2、通过课题的理论设计和在计算机中实验调试代码，加深对数值分析知识的理解，培养软件开发的实践技能，提高分析解决具体问题的能力；3、充分利用面向对象程序设计，用精要的语句设计程序，培养独立创新能力；4、熟练运用MATLAB7开发界面美观、清晰，简单，好学，好用的计算软件；说明书中一般包括任务的提出、方案论证、设计与计算说明、试验调试及结果的分析、结束语等内容。
要求理论依据充分，数据准确，公式推导及计算结果正确。
涉及到计算机软件：学生要独立完成一个软件或较大软件中的一个模块，要有足够的工作量；要写出相关的设计文件；能够进行计算机演示和给出运行结果。
在此软件中，我们运用龙贝格方法进行积分。
此外，只需要在编辑框中输入上下限和含x的被积函数，点击计算按钮，得到计算结果。
关键词：MATLAB 数值分析龙贝格积分绪论数值分析是计算数学的一个主要部分，它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。
数值分析一个很强的特点就是面向计算机，要根据计算机特点提供切实可行的有效算法。
因此，本次课程设计是在学习了《数值分析》和相关开发软件课程之后，为了能掌握数值分析的基本概念，结合实际的操作和设计，巩固课堂教学内容，使自己能掌握数值分析的基本概念、原理和技术，将理论与实际相结合，应用现有的数据建模工具和数值分析理论，规范、科学地完成一个小型可视化软件的设计与实现，把理论课与实验课所学内容做一综合，通过实际项目的设计、开发，培养学生独立进行计算软件的建模、在计算机中进行计算、设计、通过相关软件开发系统，并在此基础上强化学生的实践意识、提高其实际动手能力和创新能力。
数值积分算法与MATLAB实现 论
编号：审定成绩：毕业设计（论文）设计（论文）题目：数值积分算法与MATLAB实现学院名称：数理学院学生姓名：专业：数学与应用数学班级：学号：指导教师：答辩组负责人：填表时间：年月摘要在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。
数值积分就是解决此类问题的一种行之有效的方法。
积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。
本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。
本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。
除了研究这些数值积分算法的理论外，本文还将这些数值积分算法在计算机上通过MATLAB软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。
【关键词】数值积分牛顿-科特斯求积公式高精度求积公式 MATLAB软件ABSTRACTWhen the solution of the definite integral of some function values，because the original function is very complex and difficult to find the elementary function expression, the integral is difficult to accurately calculate, only managed to find the approximate value, and the case is small that allows to direct interface with the Newton - Leibniz formula to calculate the definite integral. Numerical integration is an effective method to solve such problems. The numerical integration is an important branch of numerical analysis; therefore, exploring the approximate calculation of the numerical integration method has obvious practical significance. This article departure from the numerical integration problem, described in detail some important numerical integration methods.This paper has introduced detail the Newton - Coates quadrature formula, and in order to improve the calculation accuracy of numerical integration formulas, More precise formulas have Romberg quadrature formulas and the Gauss - Legendre quadrature formula. In addition to the study of these numerical integration algorithm theory, the article also involve what these numerical integration algorithm be programmed by matlab software on the computer, and an example is calculated with a variety of quadrature formulas, finally analysis and comparison to various quadrature formulas calculation error.【Key words】Numerical integration Newton-Cotes quadrature formulaHigh-precisionquadrature formula Matlab software目录前言 (1)第一章牛顿-科特斯求积公式 (2)第一节数值求积公式的构造 (2)第二节复化求积公式 (9)第三节本章小结 (12)第二章高精度数值积分算法 (13)第一节梯形法的递推 (13)第二节龙贝格求积公式 (14)第三节高斯求积公式 (17)第四节高斯-勒让德求积公式 (19)第五节复化两点高斯-勒让德求积公式 (22)第六节本章小结 (23)第三章各种求积公式的MATLAB编程实现与应用 (24)第一节几个低次牛顿-科特斯求积公式的MATLAB实现 (24)第二节复化求积公式的MATLAB实现 (28)第三节龙贝格求积公式的MATLAB实现 (33)第三节高斯-勒让德求积公式的MATLAB实现 (34)第五节各种求积算法的分析比较 (36)第六节本章小结 (38)结论 (39)致谢 (40)参考文献 (41)附录 (43)一、英文原文 (43)二、英文翻译 (53)前 言对于定积分()ba f x dx ⎰，在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式()()()ba I f x dx Fb F a ==-⎰可以计算定积分的值，但在很多情况下()f x 的原函数不易求出或非常复杂。
蒙特卡罗计算积分 matlab
蒙特卡罗计算积分 matlab
蒙特卡罗计算积分是一种数值计算方法，其基本思想是通过随机抽样的方式来求解积分。
在 MATLAB 中，可以通过编写代码实现蒙特卡罗积分的计算。
具体步骤如下：
1. 定义需要求解的函数。
2. 设定积分区间和随机抽样点的个数。
3. 生成随机数，计算函数值并统计在积分区间内的点数。
4. 计算积分值。
以下是 MATLAB 代码示例：
function [I, x] = monte_carlo_int(f, a, b, n)
% f: 被积函数
% a, b: 积分区间
% n: 随机抽样点的个数
x = a + (b - a) * rand(n, 1); % 生成随机数
fx = f(x); % 计算函数值
I = (b - a) * sum(fx) / n; % 计算积分值
end
例如，要求解函数 f(x) = x^2 在区间 [0,1] 上的积分，可以使用以下代码：
f = @(x) x.^2;
a = 0;
b = 1;
n = 100000;
[I, x] = monte_carlo_int(f, a, b, n);
disp(['积分值为：', num2str(I)]);
运行结果可能会有所不同，因为随机抽样的结果是随机的。
但通常情况下，随着随机抽样点的增加，计算结果会越来越接近真实的积分值。
Romberg求积分公式
f=0.1
v0=3
w0=2
t =
0.3398
v =
2.6667
三、高程
1、题中已给出4*4=16个数据，分别对应16个坐标位置上的高程，现只需采用插值的方法，向其中填补数值，便可拟合对应的曲面，考虑到找到适合曲面的二元函数比较复杂，并且插值之后的数据量够大（10000个），具有一定的代表性，因此在求解丘陵最高点及其高程的时候，可以将所有数据进行比较，取其最大值所对应的x，y值作为最高点。具体程序如下
694.7456 696.7677 698.6979 700.5364 702.2831 703.9379 705.5010
695.6694 697.6433 699.5262 701.3184 703.0196 704.6300 706.1494
696.5197 698.4463 700.2829 702.0295 703.6862 705.2530 706.7297
w=0;
if(h~=0)
fori=1:(2^(k-1)-1)
w=w+f(a+i*h);
end
t(k,1)=h*(f(a)/2+w+f(b)/2);
forl=2: k
fori=1:(k-l+1)
t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);
end
end
658.3765 661.6218 664.7536 667.7717 670.6763 673.4673 676.1448
660.8430 664.0209 667.0864 670.0395 672.8801 675.6083 678.2240
数值积分龙贝格matlab
《数值分析》课程实验报告一、实验目的1、进一步熟悉向量矩阵的运算；2、掌握龙贝格（Romberg ）算法，并能用高级程序语言MATLAB 编写实现此算法的程序；3、进而加深对龙贝格（Romberg ）算法的理解。
二、实验内容1. 使用Romberg 积分，对于计算下列⎰+4802)cos (1dx x 各近似值a.确定1,51,41,31,21,1,,,,R R R R Rb.确定5,54,43,32,2,,,R R R Rc.6,65,64,63,62,61,6,,,,,R R R R R Rd.确定10,109,98,87,7,,,R R R R三、实验步骤1. 编写程序龙贝格积分方法如下：n=5;a=0;b=48;h(1,1)=b-a;fa=sqrt(1+(cos(a))^2);fb=sqrt(1+(cos(b))^2);r(1,1)=h(1,1)/2*(fa+fb);disp('R11，R21，R31，R41，R51分别为');disp(r(1,1));for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);disp(r(i,1));enddisp('R22，R33，R44，R55分别为');for k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);enddisp(r(k,k));enddisp('R61，R62，R63，R64，R65，R66分别为');n=6;for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);endfor k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);endendfor i=1:ndisp(r(6,i));enddisp('R77，R88，R99，R10,10分别为');n=10;for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);endfor k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);endendfor i=7:10disp(r(i,i));end运行结果如下：R11，R21，R31，R41，R51分别为62.437457.288656.443856.263156.2188R22，R33，R44，R55分别为55.572356.201556.205656.2041R61，R62，R63，R64，R65，R66分别为58.362759.077359.268959.317559.329759.3328R77，R88，R99，R10,10分别为58.422158.470758.470558.4705四、实验小结在这次编程中我学到了很多东西，把程序写进软件中也出现了很多错误，细节问题使我们必须注意的，自己有了很多的收获，自己进一步理解和学习了Matlab软件。
龙格库塔方法及其matlab实现
资料范本本资料为word版本，可以直接编辑和打印，感谢您的下载龙格库塔方法及其matlab实现地点：__________________时间：__________________说明：本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。
MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MatLab 是英文MATrix LABoratory（矩阵实验室）的缩写。
在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
关键词：龙格-库塔 matlab 微分方程前言1.1：知识背景龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。
该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。
Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。
经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。
- 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。