
#include<math.h>
using namespace std
double fun(double x)
{
return (x*sin(x)) /*修改此处可以改变被积函数*/
}
double integral(double a, double b, int n)
{
double sum,step,result
int i
如如 sum=(fun(a)+fun(b))/2
step=(b-a)/n /*积分步长*/
for(i=1i<ni++)
sum=sum+fun(a+i*step)
裤消 result=sum*step
return result/*返回积分值*/
}
void main() {
double result
result=integral(1.0,2.0,150)/*修改此处可以改变积分上下胡橡知限和步长*/
printf("result=%f\n",result)
system("pause")
}
运行结果:
result=1.440418
main.m
g=@(t)exp(-t.^2)
x=linspace(0,5,100)
y1=zeros(1,100)
y2=zeros(1,100)
for i = 1:100
y1(i)=2/sqrt(pi)*Simpson(g,0,x(i),1000) %幸普森公式
y2(i)=2/sqrt(pi)*trapezoidal(g,0,x(i),1000)%复合梯形公式
end
plot(x,y1,'r-',x,y2,'b.')
legend('辛普森公式','复合梯形公式')
Simpson.m
function y = Simpson(f,a,b,M)
% f被积函数;a积分下限;b积分上限猜桥;M子区间个数(将x分为多少个区间穗磨猛)
h=(b-a)/(2*M)
s1=0
s2=0
for i=1:M
x=a+(2*i-1)*h
s1=s1+f(x)
end
for j=1:(M-1)
x=a+2*j*h
s2=s2+f(x)
end
y=h/3*(f(a)+2*s2+4*s1+f(b))
end
trapezoidal.m
function y = trapezoidal(f,a,b,n)
h=(b-a)/n
x=zeros(1,n+1)
F=zeros(1,n+1)
for k=0:n
x(k+1)=a+k*h
if x(k+1)==0
x(k+1)=10^(-10)
游枣 end
end
T_1=h/2*(f(x(1))+f(x(n+1)))
for i=2:n
F(i)=h*f(x(i))
end
T_2 = sum(F)
y =T_1+T_2
end
欢迎分享,转载请注明来源:内存溢出
微信扫一扫
支付宝扫一扫
评论列表(0条)