2013-04-14 9 views
4

) 숫자 방법 클래스의 경우 Simpson의 합성 규칙을 사용하여 명확한 적분을 평가하는 프로그램을 작성해야합니다. 나는 이미 이것 (아래 참조)을 가지고 있지만, 나의 대답은 정확하지 않다. 나는 f (x) = x 인 프로그램을 0에서 1로 적분하여 결과를 0.5로 테스트한다. 나는 0.78746 ... 등을 얻는다. 나는 Scipy에서 이용할 수있는 Simpson 's 규칙이있다라는 것을 알고있다. 그러나 나는 정말로 그것을 직접 작성해야한다.Simpson의 규칙 (Python

두 개의 루프에 문제가 있다고 생각됩니다. 나는 "범위 (1, n, 2)에", "범위 (2, n-1, 2)에"전에 시도했는데, 이것은 0.41668333 ... 등의 결과를 줬다. 나는 또한 시도했다. "x + = h"이고 "x + = i * h"시도했습니다. 첫 번째는 0.3954, 두 번째 옵션은 7.9218입니다.

# Write a program to evaluate a definite integral using Simpson's rule with 
# n subdivisions 

from math import * 
from pylab import * 

def simpson(f, a, b, n): 
    h=(b-a)/n 
    k=0.0 
    x=a 
    for i in range(1,n/2): 
     x += 2*h 
     k += 4*f(x) 
    for i in range(2,(n/2)-1): 
     x += 2*h 
     k += 2*f(x) 
    return (h/3)*(f(a)+f(b)+k) 

def function(x): return x 

print simpson(function, 0.0, 1.0, 100) 
+2

파이썬 2 파이썬 (3) 내의 2ln (X/5)를 통합 할 것인가? – Makoto

+0

이것은 절단 오류로 인한 것일 수 있습니다. 어느 시점에서 불확실한 값을 확인 했습니까? – StoryTeller

+0

http://en.wikipedia.org/wiki/Simpsons_rule을 보셨습니까 python2에 나열된 알고리즘이 있습니다 – epsilonhalbe

답변

4

당신은 아마 또한, 시작 조건 및 반복 횟수가 꺼져, 두 번째 루프 전에 x를 초기화하는 것을 잊지. 올바른 방법은 다음과 같습니다.

def simpson(f, a, b, n): 
    h=(b-a)/n 
    k=0.0 
    x=a + h 
    for i in range(1,n/2 + 1): 
     k += 4*f(x) 
     x += 2*h 

    x = a + 2*h 
    for i in range(1,n/2): 
     k += 2*f(x) 
     x += 2*h 
    return (h/3)*(f(a)+f(b)+k) 

실수는 루프 불변의 개념과 관련이 있습니다. 세부 사항을 너무 많이 이해하지 못하는 것은 일반적으로주기가 끝나고 사이클의 끝에서 진행되는 사이클을 이해하고 디버그하는 것이 더 쉽습니다. 여기서는 x += 2 * h 행을 끝까지 이동 했으므로 합계가 시작되는 위치를 쉽게 확인할 수있었습니다. . 당신의 구현에서 루프의 첫 번째 라인으로 2 * h을 추가하는 첫번째 루프에만 이상한 x = a - h을 할당 할 필요가 있습니다.

+0

감사합니다. 그러나, 나는 수학적 실수와 더 많은 관계가 있다고 생각한다. (각 루프는 시작 값으로 각각 다른 x, 각각 x1과 x2가 필요하다.) 나는 몇 분 전에이 솔루션을 직접 사용했기 때문에 계산을했다. x는 k를 계산하는 선 위에 있습니다. 그러면 0.4997333의 해답을 얻을 수 있고 0.4802666의 해답을 얻을 수 있습니다.당신이 말하는 것을 뒤늦게 보았을 때 루프의 끝 부분에 x를 계산하는 선을 두는 것이 더 합리적이기 때문에 이것은 혼란 스럽습니다. – Geraldine

+0

또한 파이썬 질문보다 더 많은 수학 문제가 있습니다. 제 책에 따르면 Simpson 's는 사다리꼴 규칙보다 정확해야합니다. 나는 또한 사다리꼴 규칙에 대한 프로그램을 썼다. n = 100 일 때, 사다리꼴 규칙은 정확한 답을 주지만 (Simpson 's의 구현과 함께, 그것은 0.4802666입니다.) 왜 어떤 아이디어? 내 코드인가? – Geraldine

+0

자, 작동중인 위키 피 디아 코드와 광산을 가지고 있습니다. 둘 다 n = 100에 대해 0.5를줍니다. 코드가 다른 것이라면 Simpson의 규칙이 아닙니다. 잘못되었습니다. 비교할 가치가 없습니다. – unkulunkulu

1

이 코드를 작성하기 위해서해야 할 일은 bounds 함수에 a와 b에 대한 변수를 추가하고 x 변수를 사용하는 f (x)에 함수를 추가하는 것입니다. 원하는 경우 함수와 범위를 simpsonsRule 함수에 직접 구현할 수도 있습니다 ... 또한 이들은 프로그램 자체가 아닌 프로그램에 임 플린 팅되는 함수입니다.

def simpsonsRule(n): 

    """ 
    simpsonsRule: (int) -> float 
    Parameters: 
     n: integer representing the number of segments being used to 
      approximate the integral 
    Pre conditions: 
     Function bounds() declared that returns lower and upper bounds of integral. 
     Function f(x) declared that returns the evaluated equation at point x. 
     Parameters passed. 
    Post conditions: 
     Returns float equal to the approximate integral of f(x) from a to b 
     using Simpson's rule. 
    Description: 
     Returns the approximation of an integral. Works as of python 3.3.2 
     REQUIRES NO MODULES to be imported, especially not non standard ones. 
     -Code by TechnicalFox 
    """ 

    a,b = bounds() 
    sum = float() 
    sum += f(a) #evaluating first point 
    sum += f(b) #evaluating last point 
    width=(b-a)/(2*n) #width of segments 
    oddSum = float() 
    evenSum = float() 
    for i in range(1,n): #evaluating all odd values of n (not first and last) 
     oddSum += f(2*width*i+a) 
    sum += oddSum * 2 
    for i in range(1,n+1): #evaluating all even values of n (not first and last) 
     evenSum += f(width*(-1+2*i)+a) 
    sum += evenSum * 4 
    return sum * width/3 

def bounds(): 
    """ 
    Description: 
     Function that returns both the upper and lower bounds of an integral. 
    """ 
    a = #>>>INTEGER REPRESENTING LOWER BOUND OF INTEGRAL<<< 
    b = #>>>INTEGER REPRESENTING UPPER BOUND OF INTEGRAL<<< 
    return a,b 

def f(x): 
    """ 
    Description: 
     Function that takes an x value and returns the equation being evaluated, 
     with said x value. 
    """ 
    return #>>>EQUATION USING VARIABLE X<<< 
0

이 프로그램을 사용하면 Simpson의 1/3 규칙을 사용하여 한정 분모를 계산할 수 있습니다. 변수 패널의 값을 늘리면 정확도를 높일 수 있습니다. 예를 들어

import numpy as np 

def integration(integrand,lower,upper,*args):  
    panels = 100000 
    limits = [lower, upper] 
    h = (limits[1] - limits[0])/(2 * panels) 
    n = (2 * panels) + 1 
    x = np.linspace(limits[0],limits[1],n) 
    y = integrand(x,*args) 
    #Simpson 1/3 
    I = 0 
    start = -2 
    for looper in range(0,panels): 
     start += 2 
     counter = 0 
     for looper in range(start, start+3): 
      counter += 1 
      if (counter ==1 or counter == 3): 
       I += ((h/3) * y[looper]) 
      else: 
       I += ((h/3) * 4 * y[looper]) 
    return I 

:

def f(x,a,b): 
    return a * np.log(x/b) 
I = integration(f,3,4,2,5) 
print(I) 

간격이 3, 4

관련 문제