复化梯形公式和复化Simpson公式

数值计算方法上机题目3

一、计算定积分的近似值: 要求:

(1)若用复化梯形公式和复化Simpson公式计算,要求误差限???10?7,分别利用他们的余项估计对每种算法做出步长的事前估计; (2)分别利用复化梯形公式和复化Simpson公式计算定积分; (3)将计算结果与精确解比较,并比较两种算法的计算量。 1.复化梯形公式 程序:

程序1(求f(x)的n阶导数:

symsx

f=x*exp(x)%定义函数f(x) n=input('输入所求导数阶数:') f2=diff(f,x,n)%求f(x)的n阶导数

12结果1 输入n=2 f2=

2*exp(x)+x*exp(x) 程序2:

clc clear

symsx%定义自变量x

f=inline('x*exp(x)','x')%定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可 f2=inline('(2*exp(x)+x*exp(x))','x')%定义f(x)的二阶导数,输入程序1里求出的f2即可。 f3='-(2*exp(x)+x*exp(x))'%因fminbnd()函数求的是表达式的最小值,且要求表达式带引号,故取负号,以便求最大值 e=5*10^(-8)%精度要求值 a=1%积分下限

b=2%积分上限

x1=fminbnd(f3,1,2)%求负的二阶导数的最小值点,也就是求二阶导数的最大值点对应的x值 forn=2:1000000%求等分数n

Rn=-(b-a)/12*((b-a)/n)^2*f2(x1)%计算余项 ifabs(Rn)

h=(b-a)/n%求h Tn1=0

fork=1:n-1%求连加和 xk=a+k*h

Tn1=Tn1+f(xk) end

Tn=h/2*((f(a)+2*Tn1+f(b))) z=exp(2)

R=Tn-z%求已知值与计算值的差

fprintf('用复化梯形算法计算的结果Tn=') disp(Tn)

fprintf('等分数n=') disp(n)%输出等分数

fprintf('已知值与计算值的误差R=') disp(R)

输出结果显示:

用复化梯形算法计算的结果Tn= 等分数n=7019

已知值与计算值的误差R= 公式 程序:

程序1:(求f(x)的n阶导数):

symsx

f=x*exp(x)%定义函数f(x) n=input('输入所求导数阶数:') f2=diff(f,x,n)%求f(x)的n阶导数

结果1 输入n=4

f2=

4*exp(x)+x*exp(x) 程序2:

clc clear

symsx%定义自变量x

f=inline('x*exp(x)','x')%定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可 f2=inline('(4*exp(x)+x*exp(x))','x')%定义f(x)的四阶导数,输入程序1里求出的f2即可

f3='-(4*exp(x)+x*exp(x))'%因fminbnd()函数求的是表达式的最小值,且要求表达式带引号,故取负号,一边求最大值 e=5*10^(-8)%精度要求值 a=1%积分下限 b=2%积分上限

x1=fminbnd(f3,1,2)%求负的四阶导数的最小值点,也就是求四阶导数的最大值点对应的x值 forn=2:1000000%求等分数n

Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1)%计算余项 ifabs(Rn)

h=(b-a)/n%求h Sn1=0 Sn2=0

fork=0:n-1%求两组连加和 xk=a+k*h xk1=xk+h/2

Sn1=Sn1+f(xk1) Sn2=Sn2+f(xk) end

Sn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b))%因Sn2多加了k=0时的值,故减去f(a) z=exp(2)

R=Sn-z%求已知值与计算值的差

fprintf('用Simpson公式计算的结果Sn=') disp(Sn)

fprintf('等分数n=') disp(n)

fprintf('已知值与计算值的误差R=') disp(R)

输出结果显示:

用Simpson公式计算的结果Sn=

等分数n=24

已知值与计算值的误差R=

用复化梯形公式计算的结果为:,与精确解的误差为:。等分数n=7019 用复化Simpson公式计算的结果为:,与精确解的误差为:。等分数n=24 3、柯斯特公式求积分: 程序代码:

(1)function[y,Ck,Ak]=NewtonCotes(fun,a,b,n) ifnargin==1 [mm,nn]=size(fun); ifmm>=8

error('为了保证NewtonCotes积分的稳定性,最多只能有9个等距节点!') elseifnn~=2

error('fun构成应为:第一列为x,第二列为y,并且个数为小于10的等距节点!') end

xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1;

elseifnargin==4 xk=linspace(a,b,n+1); ifisa(fun,'function_handle')

fx=fun(xk); else

error('fun积分函数的句柄,且必须能够接受矢量输入!') end else

error('输入参数错误,请参考函数帮助!') end

Ck=cotescoeff(n); Ak=(b-a)*Ck; y=Ak*fx';

(2)functionCk=cotescoeff(n) fori=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

(3)functionf=intfun(t,n,k) f=1;

fori=[0:k-1,k+1:n] f=f.*(t-i); end

代码解释: function[y,Ck,Ak]=NewtonCotes(fun,a,b,n)

联系客服:779662525#qq.com(#替换为@) 苏ICP备20003344号-4