matlab实现复化Newton-Cotes公式求积分的程序应用和代码
牛顿-柯特斯公式matlab

牛顿-柯特斯公式matlab首先,让我们来了解一下数值积分的基本概念。
数值积分是通过求取一个函数在给定区间上的近似面积来计算函数的定积分。
一种常见的数值积分方法是使用插值多项式来近似函数,并在给定区间上对该多项式进行积分。
牛顿插值多项式是由一组不同的x值和对应的函数值构成的。
该多项式通过这些点来逼近函数,并可以用于在任意点上计算函数的近似值。
牛顿插值多项式的形式如下:P(x)=f[x₀]+f[x₀,x₁](x-x₀)+f[x₀,x₁,x₂](x-x₀)(x-x₁)+...其中,f[x₀]表示函数在x₀上的值,f[x₀,x₁]表示函数在x₀和x₁上的差商。
柯特斯系数用于计算牛顿插值多项式在给定区间上的积分。
公式如下:C₀=1C₁=h/2C₂=h²/6C₃=h³/12C₄=h⁴/20其中,h表示区间的宽度。
在MATLAB中,可以使用以下代码来实现牛顿-柯特斯公式:function result = newton_cotes(f, a, b, n)h=(b-a)/n;x=a:h:b;fx = f(x);coefficient = zeros(n+1, 1);coefficient(1) = 1;for i = 2:n+1coefficient(i) = coefficient(i-1) * (h^(i-1)) / factorial(i-1);endresult = sum(fx .* coefficient);end```在上面的代码中,`f`表示要积分的函数,`a`和`b`表示积分区间的起始点和结束点,`n`表示节点的数量。
首先,我们计算出节点的横坐标和对应的函数值。
然后,根据柯特斯系数的公式计算系数。
最后,将函数值与系数相乘,并求和,从而得到近似的积分值。
例如,我们要计算函数f(x) = sin(x)在区间[0, π/2]上的积分值,可以使用以下代码:a=0;b = pi/2;n=4;result = newton_cotes(f, a, b, n);disp(result);```运行该代码,将输出函数f(x)在区间[0,π/2]上的近似积分值。
实验一MATLAB计算复积分和留数

实验一MATLAB计算复积分和留数复积分是对于二重积分或多重积分的进一步延伸,通过引入复数域,可以将积分的范围从实数域扩展到复数域。
留数理论是对于复变函数的一种重要工具,它用于计算闭合曲线内的奇点的积分值。
在MATLAB中,可以使用复积分和留数函数来计算复积分和留数。
首先,我们介绍复积分的计算方法。
在MATLAB中,可以使用int函数来计算复积分。
int函数的一般形式为:z = int(fun, a, b)其中,fun是对应的被积函数,a和b分别是积分的下限和上限。
下面我们以计算复数域内的复积分为例进行说明。
假设我们要计算下面的复积分:∫(1/z)dz其中,积分路径为从复平面上的点1到点-1、在MATLAB中,可以使用int函数完成计算:syms zfun = 1/z;a=1;b=-1;res = int(fun, a, b)上述代码中,我们首先定义了被积函数fun,并指定了积分的下限a和上限b。
然后使用int函数进行计算,并将结果保存在res变量中。
代码运行后,可以得到res的值,即复积分的结果。
接下来,我们介绍留数的计算方法。
在MATLAB中,可以使用residue函数来计算留数。
residue函数的一般形式为:[r, p, k] = residue(b, a)其中,b和a分别是分子和分母的多项式系数。
下面我们以计算函数(1/z^2)的留数为例进行说明。
假设我们要计算函数(1/z^2)在z=0处的留数。
在MATLAB中,可以使用residue函数完成计算:b=[01];a=[100];[r, p, k] = residue(b, a)上述代码中,我们首先定义了多项式分子b和分母a的系数。
然后使用residue函数进行计算,并将结果保存在r、p和k变量中。
r对应留数的向量,p对应极点的向量,k对应常数项。
代码运行后,可以得到r、p和k的值,即留数的结果。
综上所述,我们介绍了在MATLAB中计算复积分和留数的方法和函数。
[matlab二分法程序代码]matlab牛顿迭代法程序代码
![[matlab二分法程序代码]matlab牛顿迭代法程序代码](https://img.taocdn.com/s3/m/5cd37268a0116c175f0e48ea.png)
[matlab二分法程序代码]matlab牛顿迭代法程序代码篇一: matlab牛顿迭代法程序代码牛顿迭代法主程序:function?[k,x,wuca,yx]?=?newtonk=1;yx1=fun;yx2=fun1;x1=x0-yx1/yx2;while?abs>tolx0=x1;yx1=fun;yx2=fun1;k=k+1;x1=x1-yx1/yx2;endk;x=x1;wuca=abs/2;yx=fun;end分程序1:function?y1=funy1=sqrt-tan;end分程序2:function????y2=fun1%函数fun的导数y2=x/)-1/) );end结果:[k,x,wuca,yx]?=?newtonk?=8x?=0.9415wuca?=4.5712e-08yx?=-3.1530e-14[k,x,wuca,yx]?=?newtonk?=243x?=NaNwuca?=NaNyx?=NaN篇二: 二十个JA V A程序代码1百分制分数到等级分数package pm;public class SwitchTest {//编写程序,实现从百分制分数到等级分数的转换////>=90 A// 80~89 B// 70~79 C// 60~69 D// public static void main {int s=87;switch{case 10 :System.out.println;break; case 9 :System.out.println;break; case 8 :System.out.println;break;case 7 :System.out.println;break;case 6 :System.out.println;break; default :System.out.println;break; }}}2成法口诀阵形package pm;public class SwitchTest{public static void main{for{for{System.out.print+”\t”); }System.out.println;}}}3华氏和摄氏的转换法package pm;import java.util.Scanner;public class SwitchTest {public static void main {Scanner sc=new Scanner;while {System.out.println;String s = sc.next.trim;if ) {//做摄氏向华摄的转换System.out.println;double db = sc.nextDouble; double db2 = + 32;System.out.println;} else if ) {//做华摄向摄氏的转换System.out.println;double db = sc.nextDouble; double db2 = * 5 / 9;System.out.println + “C”); }else if){ break;}}}}package pm;import java.util.Scanner;public class SwitchTest{public static void main {Scanner sc=new Scanner;boolean flag=true;while {System.out.println; String str = sc.nextLine.trim; if || str.endsWith) {//做摄氏向华摄的转换30cString st = str.substring - 1);double db = Double.parseDouble;//[0,2)//2 double db=Double.valueOf.doubleValue; double db2 = + 32;System.out.println;} else if || str.endsWith) {//做华摄向摄氏的转换String st = str.substring - 1);double db = Double.parseDouble;//[0,2)//2 double db=Double.valueOf.doubleValue; double db2 = * 5 / 9;System.out.println + “C”); }else if){flag=false;}}}}4三个数的最大数package pm;public class SwitchTest {public static void main {int a=1,b=2,c=3,d=0;d=a>b?a:b;d=a>b?:;System.out.println;}}5简单计算器的小程序package one;import java.awt.BorderLayout;import java.awt.GridLayout;import java.awt.event.ActionEvent; import java.awt.event.ActionListener;import javax.swing.JButton;import javax.swing.JFrame;import javax.swing.JPanel;import javax.swing.JTextField;public class Jsq implements ActionListener {private JFrame frame;private JButton[] bus;private JTextField jtx;private JButton bu;private char[] strs;private String d_one = ““;private String operator;public static void main { new Jsq;}/* 利用构造进行实例化*/ public Jsq { frame = new JFrame; jtx = new JTextField; bus = new JButton[16]; strs = “789/456*123-0.+=“.toCharArray; for { bus[i] = new JButton; bus[i].addActionListener; } bu = new JButton; bu.addActionListener; init; } /* GUI 初始化*/ public void init { JPanel jp1 = new JPanel; jp1.add; jp1.add; frame.add; } /* 事件的处理*/ public void actionPerformed { /*获取输入字符*/ String conn = arg0.getActionCommand; /*清除计算器内容*/ if ) { JPanel jp2 = new JPanel; jp2.setLayout); for { jp2.add; } frame.add; frame.pack; frame.setLocation; frame.setVisible; frame.setDefaultCloseOperation;d_one = ““; operator = ““; jtx.setText; return; } /*暂未实现该功能*/if){ return; } /*记录运算符,保存运算数字*/ if ) != -1) { if && ““.equals)) return; d_one = jtx.getText; operator = conn; jtx.setText; return; } /*计算结果*/ if ) { if && ““.equals)) return; double db = 0; if ) { db = Double.parseDouble + Double.parseDouble); jtx.setText; } if ) { db = Double.parseDouble - Double.parseDouble); jtx.setText; } if ) { db = Double.parseDouble * Double.parseDouble); jtx.setText; } if ) { db = Double.parseDouble / Double.parseDouble); jtx.setText; } d_one = db + ““; return; }//界面显示jtx.setText + conn);}}6三角形图案package pm;public class SwitchTest{public static void main{int n=5;for{for{System.out.print;}for{System.out.print;}System.out.println;}}}7输出输入的姓名package pm;import java.util.Scanner;public class SwitchTest{public static void main{String name=null;Scanner sca=new Scanner ; char firstChar; do{System.out.println; name=sca.nextLine; firstChar=name.charAt;}while);System.out.println;}}8一小时倒计时小程序package pm;import javax.swing.JFrame;import javax.swing.JLabel;import javax.swing.JPanel;public class SwitchTest {private JFrame frame;private JLabel jl1;private JLabel jl2;private JLabel jl3;/*主方法*/public static void main {new SwitchTest.getTime;}/*倒计时的主要代码块*/private void getTime{long time=1*3600;long hour =0 ;long minute =0 ;long seconds=0;while{hour=time/3600;minute=/60;seconds=time-hour*3600-minute*60; jl1.setText;jl2.setText;jl3.setText;try {Thread.sleep;} catch {e.printStackTrace;}time--;}}/*构造实现界面的开发GUI */ public SwitchTest{frame = new JFrame;jl1 = new JLabel;jl2 = new JLabel;jl3 = new JLabel;init;}/*组件的装配*/private void init{JPanel jp=new JPanel;jp.add;jp.add;jp.add;frame.add;frame.setVisible;frame.setLocation;frame.setSize;frame.setDefaultCloseOperation; } }9棋盘图案public class Sjx{public static void main{int SIZE=19;for{if{System.out.print;//两个空格}else{System.out.print);//两个空格}}System.out.println;// System.out.print:);for{System.out.print;//一个空格}else{System.out.print+” “);//一个空格}for{System.out.print;//两个空格}System.out.println;}}}10数组输出唐诗package day04;public class ArrayTest {public static void main{char[][] arr=new char[4][7];String s=“朝辞白帝彩云间千里江陵一日还两岸猿声啼不住轻舟已过万重山”; for{for{arr[i][j]=s.charAt;}for{for{System.out.print;}System.out.println;}}}11找出满足条件的最小数package day02;public class Fangk{public static void main{// for{// int q=i/1000;// int b=i/100%10;// int s=i/10%10;// int g=i%10;// if{ // System.out.println; // break; // }// }loop1: for{loop2: for{if{continue loop2;}for{for{if{ System.out.println); break loop1;}}}}}}} Min Number12判断一个数是否是素数package day02;public class Fangk{ public static void main{ int num=14;boolean flag=true;for{flag=false;break;}}if{System.out.println; }else{System.out.println; }}}////////////////////////////////////////////////////////////////////// package day04;import java.util.Scanner;public class A1{public static void main{int n;Scanner sca=new Scanner;System.out.println; n=sca.nextInt;if){System.out.println; }else{System.out.println;}public static boolean isPrimeNumber{ for{if{return false;}}return true;}}13一个数倒序排列package day02;public class Daoxu{public static void main{ int olddata=3758;int newdata=0;while{for{newdata=newdata*10+olddata%10; olddata=olddata/10; }System.out.println;}}}14将一个整数以二进制输出package day04;import java.util.Scanner; public class ArrayTest { public static void main{ int n; Scanner s=new Scanner; System.out.println; n=s.nextInt; for{if)!=0){System.out.print;}else{System.out.print;}if%8==0){System.out.print;}}}}15矩形图案package day02;public class Fangk {public static void main{ int m=5,n=6; for{System.out.print;}System.out.println;for{System.out.print;for{System.out.print; }System.out.print;System.out.println;}for{System.out.print;}}}16猜数字package day02;import java.util.Scanner;public class Csz {public static void main {Scanner s = new Scanner; int num = * 1000); int m=0; for{System.out.println; m=s.nextInt;if{System.out.println;}else if{System.out.println;}else{System.out.println; break;}if{System.out.println; }}if{System.out.println; }}}17.HotelManagerpackage hotel;import java.util.Scanner;public class HotelManager {private static String[][] rooms;// 表示房间public static void main {rooms = new String[10][12];String comm;// 表示用户输入的命令for {for {rooms[i][j] = “EMPTY”;}}//while {System.out.println;Scanner sca = new Scanner; System.gc;comm = sca.next;if ) {search;} else if ) {int roomNo = sca.nextInt;String name = sca.next;in;} else if ) {int roomNo = sca.nextInt;out;} else if ) {System.out.println;break;} else {System.out.println; }}}private static void out {if-1][-1])){ System.out.println;} return; } rooms[-1][-1]=“EMPTY”; System.out.println; private static void in { if-1][-1])){ System.out.println;return;}rooms[-1][-1]=name;System.out.println;}private static void search {for {//打印房间号for {if {System.out.print + “ “); } else {System.out.print + “ “); }}//打印房间状态System.out.println;for {System.out.print;}System.out.println;}}}18.StudentManagerpackage day05.student_manager;import java.util.Scanner;public class StudentManager {static int[][] scores=new int[6][5];static String[] students={“zhangsan”,”lisi”,”wangwu”,”zhaoliu”,”qianqi”,”liuba”}; static String[] courses={“corejava”,”jdbc”,”servlet”,”jsp”,”ejb”};public static void main {for{for{scores[i][j]=*100);}}Scanner s=new Scanner; String comm;while{System.out.println; comm=s.next;if){String para=s.next; avg;}else if){String course=s.next; sort;}else if){String student=s.next; String course=s.next; get;}else if){break;}else{System.out.println; }}}//main end!public static void avg{int sIndex=-1;//int cIndex=-1; for{ if){ sIndex=i; } } if{ for{ if){ cIndex=i; } } } if{ System.out.println; return; } double avg=0.0; if{ for{ avg+=scores[sIndex][i]; } avg/=scores[sIndex].length; System.out.println; }else{ for{ avg+=scores[i][cIndex]; } avg/=scores.length; System.out.println; } } public static void sort{ int[] courseScore=new int[scores.length]; if){//如果求总分的排名//求出每个学生的总分,将成绩存放在courseScore数组中for{ int studentSum=0; for{ studentSum+=scores[i][j]; }courseScore[i]=studentSum; } }else{//如果不是求总分排名int cIndex=-1; for{//找到这门课程的下标if){ cIndex=i; } } if{//如果是一门有效的课程//把scores数组中这一列的值放到courseScore数组中!for{ courseScore[i]=scores[i][cIndex]; } }else{//如果不是一门有效的课程System.out.println; return; } } String[] studentCopy=new String[students.length]; System.arraycopy; for{ for{ if{ int temp=courseScore[i]; courseScore[i]=courseScore[j]; courseScore[j]=temp; String stemp=studentCopy[i];studentCopy[i]=studentCopy[j]; studentCopy[j]=stemp; } } } int order=1; System.out.println; for{ if{ order--; }else{ order=i+1;} System.out.print;System.out.print;System.out.println;order++;}}public static void get{int sIndex=-1;int cIndex=-1;for{if){sIndex=i;}}if{System.out.println;return;}if){//如果求总分int studentSum=0;for{studentSum+=scores[sIndex][j];}System.out.println; return;}for{if){cIndex=i;}}if{System.out.println;return;}System.out.println;}}19.Fivepackage hotel;import java.util.Scanner;/*** 首先在程序第一次运行的时候,构建出棋盘,切以后* 不能再从新构建,知道结束,所以将其放到静态代码块中。
matlab用牛顿柯特斯公式计算积分

牛顿-柯特斯公式是数值分析中常用的积分计算方法,特别适用于对函数在一定区间上的定积分进行近似计算。
在MATLAB中,我们可以利用牛顿-柯特斯公式来进行积分计算,从而获得函数在给定区间上的近似积分值。
让我们来理解一下牛顿-柯特斯公式的基本原理。
牛顿-柯特斯公式的核心思想是利用一系列的节点和相应的权重来逼近被积函数,从而得到积分的近似值。
在MATLAB中,我们可以通过内置的函数或自定义函数来实现牛顿-柯特斯公式的计算。
在使用MATLAB计算积分时,我们首先需要确定被积函数的表达式以及积分的区间。
我们可以选择合适的牛顿-柯特斯公式来进行计算。
MATLAB提供了多种内置的积分计算函数,例如quad和integral等,它们可以方便地实现对定积分的计算。
除了使用内置函数,我们还可以编写自定义的牛顿-柯特斯公式计算程序。
这样可以更灵活地控制节点和权重的选择,从而得到更精确的积分近似值。
编写自定义的牛顿-柯特斯公式计算程序可以加深对该方法的理解,并且在特定问题上可能获得更好的计算结果。
在实际应用中,牛顿-柯特斯公式可以广泛用于工程、科学和数学等领域。
在信号处理中,我们可以利用牛顿-柯特斯公式对信号的频谱进行积分近似计算;在物理学中,我们可以利用牛顿-柯特斯公式对连续介质的密度分布进行积分近似计算。
牛顿-柯特斯公式的灵活性和高效性使得它成为了数值分析中不可或缺的工具。
回顾本文,我们首先介绍了牛顿-柯特斯公式的基本原理,然后讨论了在MATLAB中如何利用内置函数或自定义函数来实现积分的计算。
我们还探讨了牛顿-柯特斯公式在实际应用中的广泛性和重要性。
通过本文的阐述,我们希望读者能够更深入地理解牛顿-柯特斯公式的计算方法,并且能够灵活运用于自己的问题当中。
在个人观点和理解方面,我认为牛顿-柯特斯公式作为一种数值积分计算方法,具有较高的精度和灵活性,能够有效地解决实际问题中的积分计算需求。
在MATLAB中,利用牛顿-柯特斯公式进行积分计算不仅简单方便,而且还能获得较为准确的结果。
matlab实现数值分析插值及积分

数值分析学院:计算机专业:计算机科学与技术班级:xxx学号:xxx姓名:xxx指导教师:xxx数值分析摘要:数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科,是数学的一个分支,它以数字计算机求解数学问题的理论和方法为研究对象。
在实际生产实践中,常常将实际问题转化为数学模型来解决,这个过程就是数学建模。
学习数值分析这门课程可以让我们学到很多的数学建模方法。
分别运用matlab数学软件编程来解决插值问题和数值积分问题。
题目中的要求是计算差值和积分,对于问题一,可以分别利用朗格朗日插值公式,牛顿插值公式,埃特金逐次线性插值公式来进行编程求解,具体matlab代码见正文。
编程求解出来的结果为:P4x=x4+1。
其中Aitken插值计算的结果图如下:对于问题二,可以分别利用复化梯形公式,复化的辛卜生公式,复化的柯特斯公式编写程序来进行求解,具体matlab代码见正文。
编程求解出来的结果为:0.6932其中复化梯形公式计算的结果图如下:问题重述问题一:已知列表函数表格分别用拉格朗日,牛顿,埃特金插值方法计算P4x 。
问题二:用复化的梯形公式,复化的辛卜生公式,复化的柯特斯公式计算积分121xdx ,使精度小于5×10-5。
问题解决问题一:插值方法对于问题一,用三种差值方法:拉格朗日,牛顿,埃特金差值方法来解决。
一、拉格朗日插值法:拉格朗日插值多项式如下:首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数,对任一点i x 所对应的插值基函数)(x l i ,由于在所有),,1,1,,1,0(n i i j x j +-=取零值,因此)(x l i 有因子)())(()(110n i i x x x x x x x x ----+- 。
又因)(x l i 是一个次数不超过n 的多项式,所以只可能相差一个常数因子,固)(x l i 可表示成:)())(()()(110n i i i x x x x x x x x A x l ----=+-利用1)(=i i x l 得:)())(()(1110n i i i i i i x x x x x x x x A ----=+-于是),,2,1,0()())(()()())(()()(110110n i x x x x x x x x x x x x x x x x x l n i i i i i i n i i i =--------=+-+-因此满足i i n y x L =)( ),2,1,0(n i =的插值多项式可表示为:∑==nj j j n x l y x L 0)()(从而n 次拉格朗日插值多项式为:),,2,1,0()()(0n i x l y x L nj i j j i n ==∑=matlab 编程:编程思想:主要从上述朗格朗日公式入手:依靠循环,运用poly ()函数和conv ()函数表示拉格朗日公式,其中的poly (i )函数表示以i 作为根的多项式的系数,例如poly (1)表示x-1的系数,输出为1 -1,而poly (poly (1))表示(x-1)*(x-1)=x^2-2*x+1的系数,输出为1 -2 1;而conv ()表示多项式系数乘积的结果,例如conv (poly (1),poly (1))输出为1 -2 1;所以程序最后结果为x^n+x^n-1+……+x^2+x+1(n 的值据结果的长度为准)的对应系数。
复合辛普森公式matlab例题

文章标题:深度剖析复合辛普森公式在Matlab中的应用1.引言复合辛普森公式是数值分析中常用的积分逼近方法,在工程和科学领域都有广泛的应用。
本文将深入探讨复合辛普森公式的原理和在Matlab中的具体例题应用,以帮助读者全面理解该方法的实际操作和应用场景。
2.复合辛普森公式简介复合辛普森公式是一种数值积分方法,通过将积分区间分割成若干个小区间,然后在每个小区间上使用辛普森公式进行积分逼近,最终得到整个积分区间上的近似值。
其公式表达为:\[ S_n(h) = \frac{h}{3}[f(x_0) + 4\sum_{i=1}^{n/2}f(x_{2i-1}) +2\sum_{i=1}^{n/2-1}f(x_{2i}) + f(x_n)] \]其中,\(h\)为步长,\(n\)为分割的小区间数。
3. Matlab实例应用假设要对函数\(f(x) = x^3 + 2x^2 + 4x + 1\)在区间\([a, b]\)上进行积分逼近,可以通过Matlab编程实现复合辛普森公式的应用。
需要确定积分区间的上下限,然后计算步长\(h\),接着编写Matlab代码进行求解。
```matlabfunction result = simpson(f, a, b, n)h = (b - a) / n;x = a:h:b;y = f(x);result = h/3 * (y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end));end% 调用simpson函数进行积分逼近f = @(x) x.^3 + 2*x.^2 + 4*x + 1;a = 1;b = 2; n = 4;result = simpson(f, a, b, n);disp(result);```在以上Matlab代码中,首先定义了一个名为simpson的函数,用于计算复合辛普森公式的近似积分值。
Newton-Cotes梯形公式数值积分及其MATLAB范例

n I error 3 0. 555143410 - 0. 404428908 9 0. 926771814 - 0. 032800504 33 0. 957485472 - 0. 002086846 65 0. 859050139 - 0. 000522179 从输出结果可以看出, 误差 error 随着 n 的增加而减少。 更多范例限于篇幅不再赘述。 参考文献: [ 1] 曾喆昭,王耀南一种高精度数值积分方法 [J] . 湖南大 学学报,2007 年 1 月 [ 2] Gerala Reckttenwald ( 美) . 数值方法和 Mtalab 实现与应用 [ M] . 北京: 机械工业出版社,2005. [ 3] 黄明游. 数值计算方法 [ M] . 北京: 科学出版社,2005. 基金项目: 海南省自然科学基金资助项目 ( 编号: 112008 )
∫ xe
0
5
-x
dx ,将结果列出来:
图2
梯形公式应用与被 n 个点划分为 n - 1 个小段的积分区间
3 、任意 f ( x) 的梯形公式的 Matlab 实现和范例 程序清单中的 trapezoid 函数使用复合梯形公式对用户自定 义函数进行数值积分 。函数 trapezoid 的调用方法为: I = trapezoid ( fun,a,b,npanel) 其中 fun 是求被积函数的 m 文件的名字, a 和 b 是积分的 上界和下界,npanel 是积分中所使用小段的数目 。 数值积分的 精度受 npanel 的影响很大,所以一般在应用 trapezoid 时要尝试 多个 npanel 值来调用函数以使积分值落在期望的容差内 。 程序清单函数 trapezoid 使用梯形公式计算定积分 function I = trapezoid ( fun,a,b,npanel) % trapezoid Composite trapezoid rule % Synopsis: I = trapezoid ( fun,a,b,npanel)
用MATLAB计算某些区域上的二重积分

用MATLAB计算某些区域上的二重积分摘要:本文研究某些区域上二重积分数值积分公式的构造及用数学软件MATLAB 实现所构造的数值积分的计算,通过MATLAB 用所得公式计算某些典型的二重积分。
主要工作包括:将定积分数值计算的几个公式推广到二重积分,编制MATLAB 程序,最后通过具体的数值算例进行精度比较,从中选出精度高的二重积分的计算公式,并用公式计算一些典型的二重积分。
关键词:二重积分;数值计算;插值多项式;求积公式 1.引言二重积分的计算在科学计算中起着重要的作用,关于矩形区域上的二重积分的计算一般都是化重积分为累次积分,然后借助定积分已有的数值积分计算公式推导出,MATLAB 已经有这些计算公式的相应的函数,但是往往我们建模得到的二重积分的积分区域都不是矩形区域,对于一般的非矩形区域的二重积分,直接用MATLAB 是无法计算的。
又当被积函数比较复杂,无法用初等函数表示或求其原函数很困难时,就只能求积分的数值解。
若),(y x f 在D 上连续,二重积分y x y x f Dd d ),(⎰⎰存在且为一确定的常数,这个数值与),(y x f 的结构、D 的几何形状有关,二重积分计算的基本途径是在一定条件下化为二次积分,本文研究的某些区域的二重积分,要求二重积分在该区域上能化为二次积分。
二重积分的存在性[1]:),(y x f 在闭区域D 上连续,则(,)Df x y dxdy ⎰⎰必存在。
定理[1]:若),(y x f 在闭区域D {;()()}≤≤≤≤a x b c x y d x 上连续,且()c x 、()d x 在[,]a b 上连续,则:()()(,)(,)b d x ac x DI f x y dxdy f x y dxdy ==⎰⎰⎰⎰上式右端是一个先对x 后对y 的二次积分:先把),(y x f 看作x 的函数,在区间[(),()]c x d x 上对x 计算定积分(这时y 看作常数),把得到的结果(是y 的函数)再在[,]a b 上对y 计算定积分即为二重积分。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
执行函数为mymulNewtonCotes.m
1、使用方法:
Step1:在MATLAB 命令窗口输入被积函数212
0t t e dt ⎰。
输入应为:ft=@(t)t.*exp(t^2/2)。
Step2:执行函数。
输入形式为mymulNewtonCotes(ft,a,b,m,n);
其中ft —被积函数,此体重ft=@(t)t.*exp(t^2/2),已经在第一步赋值;
a —积分下限,本题中为0;
b —积分上限,本题中为1;
m —将区间[a,b]等分的子区间数量,本题可选为10;
n —采用的Newton-Cotes 公式的阶数,必须满足n<8,否则积分没法保
证稳定性。
当n=1时,即为复化梯形公式;n=2时,即为复化复化辛普森公式。
所以,分别输入mymulNewtonCotes(ft,0,1,10,1)和
mymulNewtonCotes(ft,0,1,10,2)就可以得到两种方法的积分计算结果。
2、计算结果
而根据积分运算,可得:
221112
1102222
2
0000() 1.648710.64872t t x x t t e dt e d e dx e e e ====-=-=⎰⎰⎰ 说明复化梯形和复化辛普森公式计算出的结果基本一致,与实际结果相符。
3、程序代码
function yy = mymulNewtonCotes(ft,a,b,m,n)
% 复化Newton-Cotes数值积分公式,即在每个子区间上使用Newton-Cotes公式,然后求和, % 参考的输入形式为mymulNewtonCotes(ft,0,1,10,2)
% 参数说明:
% ft——被积函数,此题中ft=@(t)t.*exp(t^2/2)
% a——积分下限
% b——积分上限
% m——将区间[a,b]等分的子区间数量
% n——采用的Newton-Cotes公式的阶数,必须满足n<8,否则积分没法保证稳定性
% (1)n=1时为复化梯形公式
% (2)n=2时为复化辛普森公式
xx = linspace(a,b,m+1);
for l = 1:m
s(l) = myNewtonCotes(ft,xx(l),xx(l+1),n);
end
yy = sum(s);
function [y,Ck,Ak] = myNewtonCotes(ft,a,b,n)
% 牛顿-科特斯数值积分公式
% Ck——科特斯系数
% Ak——求积系数
% y——牛顿-科特斯数值积分结果
xk = linspace(a,b,n+1);
for j = 1:n+1
ff(j) = ft(xk(j));
end
% 计算科特斯系数
for i=1:n+1
k=i-1;
Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n); end
% 计算求积系数
Ak=(b-a)*Ck;
% 求和算积分
y=Ak*ff';
function f=intfun(t,n,k)
% 科特斯系数中的积分表达式
f=1;
for i=[0:k-1,k+1:n]
f=f.*(t-i);
end。