# 利用辛普森 理則 (Simpson's Rule)
#計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積
from math import *
def F(x):
temp=pow(pow((1+pow(x,3)),1.0/3),2)
return temp
a=0.0
b=2.0
n=10
sum1=0.0
sum2=0.0
sn=0.0
m=int(n/2)
h=(b-a)/n
for i in range (1,2*m):
x=a+i*h
if(i%2==0):
sum2=sum2+F(x)
else:
sum1=sum1+F(x)
sn= pi *(h/3.0)*(F(a)+F(b)+2.0*sum2+4.0*sum1)
print("S{%d}={%f}\n" % (n,sn))
===============================
輸出結果======== RESTART: H:/2018-09勤益科大數值分析/數值分析/PYTHON/EX4-5.py ===========
S{10}={1.957852}
S{10}={2.394418}
S{10}={3.348841}
S{10}={3.900651}
S{10}={5.230509}
S{10}={6.048317}
S{10}={8.068276}
S{10}={9.308714}
S{10}={12.325078}
>>>
====================
辛普森理則的python另一寫法
====================
# Python code for simpson's 1 / 3 rule
import math
# Function to calculate f(x)
def func( x ):
return (math.pi* math.pow((1+math.pow(x,3)),(2/3)))
# Function for approximate integral
def simpsons_( ll, ul, n ):
# Calculating the value of h
h = ( ul - ll )/n
# List for storing value of x and f(x)
x = list()
fx = list()
# Calcuting values of x and f(x)
i = 0
while i<= n:
x.append(ll + i * h)
fx.append(func(x[i]))
i += 1
# Calculating result
res = 0
i = 0
while i<= n:
if i == 0 or i == n:
res+= fx[i]
elif i % 2 != 0:
res+= 4 * fx[i]
else:
res+= 2 * fx[i]
i+= 1
res = res * (h / 3)
return res
# Driver code
lower_limit = 0.0 # Lower limit
upper_limit = 2.0 # Upper limit
n = 10 # Number of interval
print("%.6f"% simpsons_(lower_limit, upper_limit, n))
12.325078
>>>
沒有留言:
張貼留言