复合梯形公式、复合辛普森公式 matlab | 您所在的位置:网站首页 › 复合梯形公式例题分析计算 › 复合梯形公式、复合辛普森公式 matlab |
1. 用1阶至4阶Newton-Cotes公式计算积分
程序: function I = NewtonCotes(f,a,b,type) % syms t; t=findsym(sym(f)); I=0; switch type case 1, I=((b-a)/2)*(subs(sym(f),t,a)+subs(sym(f),t,b));
case 2, I=((b-a)/6)*(subs(sym(f),t,a)+4*subs(sym(f),t,(a+b)/2)+... subs(sym(f),t,b));
case 3, I=((b-a)/8)*(subs(sym(f),t,a)+3*subs(sym(f),t,(2*a+b)/3)+... 3*subs(sym(f),t,(a+2*b)/3)+subs(sym(f),t,b));
case 4, I=((b-a)/90)*(7*subs(sym(f),t,a)+... 32*subs(sym(f),t,(3*a+b)/4)+... 12*subs(sym(f),t,(a+b)/2)+... 32*subs(sym(f),t,(a+3*b)/4)+7*subs(sym(f),t,b));
case 5, I=((b-a)/288)*(19*subs(sym(f),t,a)+... 75*subs(sym(f),t,(4*a+b)/5)+... 50*subs(sym(f),t,(3*a+2*b)/5)+... 50*subs(sym(f),t,(2*a+3*b)/5)+... 75*subs(sym(f),t,(a+4*b)/5)+19*subs(sym(f),t,b));
case 6, I=((b-a)/840)*(41*subs(sym(f),t,a)+... 216*subs(sym(f),t,(5*a+b)/6)+... 27*subs(sym(f),t,(2*a+b)/3)+... 272*subs(sym(f),t,(a+b)/2)+... 27*subs(sym(f),t,(a+2*b)/3)+... 216*subs(sym(f),t,(a+5*b)/6)+... 41*subs(sym(f),t,b));
case 7, I=((b-a)/17280)*(751*subs(sym(f),t,a)+... 3577*subs(sym(f),t,(6*a+b)/7)+... 1323*subs(sym(f),t,(5*a+2*b)/7)+... 2989*subs(sym(f),t,(4*a+3*b)/7)+... 2989*subs(sym(f),t,(3*a+4*b)/7)+... 1323*subs(sym(f),t,(2*a+5*b)/7)+... 3577*subs(sym(f),t,(a+6*b)/7)+751*subs(sym(f),t,b)); end
syms x f=exp(-x).*sin(x); a=0;b=2*pi; I = NewtonCotes(f,a,b,1)
N=1: I = 0 N=2: I = 0 N=3: I = (pi*((3*3^(1/2)*exp(-(2*pi)/3))/2 - (3*3^(1/2)*exp(-(4*pi)/3))/2))/4 N=4: I = (pi*(32*exp(-pi/2) - 32*exp(-(3*pi)/2)))/45
2. 已知,因此可以通过数值积分计算的近似值。 (1)分别取和,利用复合梯形公式和复合Simpson公式计算的近似值; 程序: function Y= CombineTraprl(f,a,b,h) %用复合梯形公式计算积分 syms t; t= findsym(sym(f)); n=(b-a)/h; I1= subs(sym(f),t,a); l=0; for k=1:n-1 xk=a+h*k; l=l+2*subs(sym(f),t,xk); end Y=(h/2)*(I1+l+subs(sym(f),t,b));
syms x f=4/(1+x^2); a=0;b=1; y= CombineTraprl(f,a,b,0.1); vpa(y,6) h=0.1: ans = 3.13993 H=0.2: ans =
1.04498 复合辛普森: function Y= CombineSimpson(f,a,b,h) %用复合辛普森公式计算积分 syms t; t= findsym(sym(f)); n=(b-a)/h; I1= subs(sym(f),t,a); l=0; for k=1:n-1 xk=a+h*k; l=l+2*subs(sym(f),t,xk); end l2=0; for k=1:n-1 xk2=a+h*(k+1)/2; l2=l2+4*subs(sym(f),t,xk2); end Y=(h/6)*(I1+l+l2+subs(sym(f),t,b));
H=0.1: ans =
3.22605 H=0.2: ans =
2.93353
(2)把区间[0,1] 等分,利用复合梯形公式和复合Simpson公式计算的近似值,若要求误差不超过,问需要把区间[0,1]划分成多少等份; function n=trap(f,a,b) syms t; t= findsym(sym(f)); I=zeros(1,500); I(1)=((b-a)/2)*(subs(sym(f),t,a)+subs(sym(f),t,b)); I(2)=((b-a)/4)*(subs(sym(f),t,a)+2*subs(sym(f),t,(b-a)/2)+subs(sym(f),t,b)); k=3; while((I(k-1)-I(k-2))>1/2*10^(-6)) l=0; for i=1:k-1 xi=a+(b-a)/k*i; l=l+2*subs(sym(f),t,xi); end I(k)=((b-a)/(2*k))*(subs(sym(f),t,a)+l+subs(sym(f),t,b)); k=k+1; end n=k-1;
syms x; f=4./(1+x.^2); a=0;b=1; n=trap(f,a,b) n =
88 复合辛普森公式: function n=Simpson(f,a,b) syms t; t= findsym(sym(f)); I=zeros(1,500); I(1)=((b-a)/6)*(subs(sym(f),t,a)+4*subs(sym(f),t,(b-a)/2)+subs(sym(f),t,b)); I(2)=((b-a)/12)*(subs(sym(f),t,a)+4*subs(sym(f),t,(b-a)/4)+4*subs(sym(f),t,3*(b-a)/4)+2*subs(sym(f),t,(b-a)/2)+subs(sym(f),t,b)); k=3; while((I(k-1)-I(k-2))>1/2*10^(-6)) l=0; m=4*subs(sum(f),t,(a+((a+b)/(2*k)))); for i=1:k-1 xi=a+(b-a)/k*i; l=l+2*subs(sym(f),t,xi); end for j=1:k-1 xj=a+(b-a)/(k*2)+(b-a)/k*j; m=m+4*subs(sym(f),t,xj); end I(k)=((b-a)/(2*k))*(subs(sym(f),t,a)+l+m+subs(sym(f),t,b)); k=k+1; end n=k-1;
n =
5 (3)选择不同的,对两种复合求积公式,试将误差描述为的函数,并比较两种方法的精度。 复合求积公式: function y=traprls(f,a,b,h) syms t; t= findsym(sym(f)); n=(b-a)/h; l=0; for k=1:n-1 xk=a+h*k; l=l+2*subs(sym(f),t,xk); end I1=(h/2)*(subs(sym(f),t,a)+l+subs(sym(f),t,b));
h=(b-a)/(n-1); n=(b-a)/h; l=0; for k=1:n-1 xk=a+h*k; l=l+2*subs(sym(f),t,xk); end I2=(h/2)*(subs(sym(f),t,a)+l+subs(sym(f),t,b));
y=I2-I1; y=abs(y); y=vpa(y,8);
syms x; f=4./(1+x.^2); a=0;b=1; h=0.01:0.05:0.5; v=zeros(1,10); for i=1:10 v(i)=traprls(f,a,b,h(i)) end v plot(h,v,'r-')
复合辛普森公式: function y=Simpsons(f,a,b,h) syms t; t= findsym(sym(f)); n=(b-a)/h; l=0; m=4*subs(sum(f),t,(a+h/2)); for k=1:n-1 xk=a++h*k; l=l+2*subs(sym(f),t,xk); end for i=1:n-1 xi=a+h/2+h*i; m=m+4*subs(sym(f),t,xi); end I1=(h/6)*(subs(sym(f),t,a)+l+m+subs(sym(f),t,b));
h=(b-a)/(n-1); n=(b-a)/h; l=0; m=4*subs(sum(f),t,(a+h/2)); for k=1:n-1 xk=a++h*k; l=l+2*subs(sym(f),t,xk); end for i=1:n-1 xi=a+h/2+h*i; m=m+4*subs(sym(f),t,xi); end I2=(h/6)*(subs(sym(f),t,a)+l+m+subs(sym(f),t,b));
y=abs(I2-I1); y=vpa(y,10);
通过图像对比可知,复合辛普森公式精度更高。 (4)是否存在某个值,当小于这个值之后,再继续减小,计算结果不再有改进?为什么? 是 复合求积公式: syms x; f=4./(1+x.^2); a=0;b=1; h=0.001:0.004:0.2; v=zeros(1,10); for i=1:50 v(i)=traprls(f,a,b,h(i)); end plot(h,v,'r-')
复合辛普森公式:
通过图像可以发现,当h |
CopyRight 2018-2019 实验室设备网 版权所有 |