%main
clc,clear;
a=0;b=1;
n=10;
trap(fun,a,b,n)
复式梯形公式
function t=trap(fun,a,b,n)
%用途:用复化梯形公式求积分
%格式:s=trap(fun,a,b,n).fun为被积函数,a,b为积分区间的左右端点
%n为区间的等分数,s返回数值积分值
syms x;
fun=@(x)(x)^(1/2)*log(x);
h=(b-a)/n;
sum=0;
for k=1:n-1
x=a+k*h;
sum=sum+feval(f,x);
end
t=h/2*(feval(f,a)+feval(f,b)+2*sum);
复式辛普森公式
function s=simp(fun,a,b,n)
%用途:用复化Simpson公式求积分
%格式:s=simp(fun,a,b,n).fun为被积函数,a,b为积分区间的左右端点
%n为区间的等分数,s返回数值积分值
syms x;
fun=@(x)(x)^(1/2)*log(x);
h=(b-a)/n;s1=0;s2=0;
for k=1:(n-1)
x=a+k*h;
s1=s1+feval(fun,x);
end
for k=0:(n-1)
x=a+h*(k+1/2);
s2=s2+feval(fun,x);
end
s=h/6*(feval(fun,a)+feval(fun,b)+2*s1+4*s2);
用变步长梯形公式计算积分,使其误差小于 10^(-4)
function [T,n]=bbctx(f,a,b,eps)
%f表示被积函数句柄;a,b表示被积区间的端点;eps表示精度
%T是用变步长梯形法求得的积分值;n表示二分区间的次数;TOL表示误差
h=b-a;fa=feval(f,a);fb=feval(f,b);T1=h*(fa+fb)/2;T2=T1/2+h*feval(f,a+h/2)/2;
n=1;
while abs(T2-T1)>=eps
h=h/2;T1=T2;S=0;x=a+h/2;
while x<b
fx=feval(f,x);S=S+fx;x=x+h;
end
T2=T1/2+S*h/2;n=n+1;
end
T=T2;
fenxi_s=int('(x)^(1/2)*log(x)',0,1);
s=vpa(fenxi_s,16);
TOL=T-s;
end
clear;clc;f=@bbctx_f;
[T,n,TOL]=bbctx(f,0,1,1e-4)
romberg.m
function [R,k,T]=romberg(fun,a,b,tol)
% 龙贝格(Romberg数值求解公式)
% inputs:
% -fun:积分函数句柄
% -a/b:积分上下限
% -tol:积分误差
% Outputs:
% -R:7阶精度Romberg积分值
% -k:迭代次数
% -T:整个迭代过程
%
% Example
% fun=@(x)(x)^(1/2)*log(x);
% [R,k,T]=romberg(fun,0,1,1e-4)
%
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a;
T=h/2*(fun(a)+fun(b));
err=1;
while err>=tol
k=k+1;
h=h/2;
tmp=0;
for i=1:n
tmp=tmp+fun(a+(2*i-1)*h);
end
T(k+1,1)=T(k)/2+h*tmp;
for j=1:k
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
end
n=n*2;
err=abs(T(k+1,k+1)-T(k,k));
end
R=T(k+1,4);
compute.m
%计算((x)^(1/2)*log(x))在0到1上面的积分
a = 0
b = 1
epsilon = 1e-4
f = @(x)(x)^(1/2)*log(x);
y = romberg(f,a,b,epsilon)
%后面是画出函数图像,注,不是积分函数图像。是被积函数图像
x = 0:0.01:1;
z = (x)^(1/2)*log(x);
plot(x,z),xlabel('x'),ylabel('y'),title(['Result = ',num2str(y)])
%自适应求积法
function I=SmartSimpson(f,a,b,eps)
if(nargin==3)
eps=1e-4;
end;
e=5*eps;
I=SubSmartSimpson(f,a,b,e);
function q=SubSmartSimpson(f,a,b,eps)
QA=IntSimpson(f,a,b,1,eps);
QLeft=IntSimpson(f,a,(a+b)/2,1,eps);
QRight=IntSimpson(f,a,(a+b)/2,1,eps);
if(abs(QLeft+QRight-QA)<=eps)
q=QA;
else
q=SubSmartSimpson(f,a,(a+b)/2,eps)+SubSmartSimpson(f,(a+b)/2,b,eps);
end
- 1
- 2
- 3
前往页