Приближение синуса и остатка

Question

Целью упражнения является оценка sin с использованием формулы ряда Маклорена.

#include <stdio.h> 
#include <math.h> 

double factorial(int n); 

int main(void) { 

    double x, p, r; 
    int i, n; 

    printf("Enter a positive double number and a non - negative integer : \n"); 
    scanf("%lf%d", &x, &n); 

    if (x <= 0) { 
     printf("Error: first argument must be a positive double.\n"); 
     return -1; 
    } 

    if (n < 0) { 
     printf("Error: second argument must be a non - negative integer.\n"); 
     return -1; 
    } 

    p = 0; 
    for (i = 0; i <= n; i++) 
    { 
     p += pow(-1, i)/factorial(2 * i + 1) * pow(x, 2 * i + 1); 
    } 

    r = fabs(sin(x) - p); 

    printf("The %d-th order Maclaurin polynomial function at x=%f is %f, with an error approximation of %f.\n", n, x, p, r); 



    getch(); 
    return 0; 
} 

double factorial(int n) 
{ 
    int i; 
    long result = 1; 
    for (i = 1; i <= n; i++) 
     result *= i; 

    return result; 
} 

Я получаю странный результат для ввода «12 16». Зачем?

Answer 1

Здесь есть три вопроса.

• Как ответил Мукит Чоудхури, долго не может содержать большие факториалы. Однако, как заметил Грегс, не должно быть никаких проблем с 16 !. Вы должны ожидать странных результатов, когда вы используете длинный, чтобы держать 21! или больше, но вам не нужно менять факториальную функцию для этого входа. Вы, вероятно, следует использовать что-то вроде следующего вместо этого, хотя:

  double factorial(int n) 
  { 
      int i; 
      double result = 1; 
      for (i = 1; i <= n; i++) 
       result *= i; 
      return result; 
  } 
  

• На входе «12 16», код утверждает, что вычисления шестнадцатого порядка Маклорена полинома, но он вычисляет 33-го порядка Маклорена многочлен вместо. Полином 16-го порядка имеет члены с точностью до -x^15/15! + 0x^16. Один из способов исправить это исправить цикл следующим образом:

  for (i = 1; i <= n; i+=2) 
  { 
      p += pow(-1, (i-1)/2)/factorial(i) * pow(x, i); 
  } 
  

  Итак, ваш код работает в проблемы с факториала, но только потому, что вы расчета дополнительных условий. Если вы вычисляете термины до -x^15/15 !, вы должны иметь возможность правильно вычислить значение полинома.

• Фактическое значение полинома Маклаурина 16-го порядка при 12 составляет -4306.756 ... что может и не быть тем, что вы ожидаете, и это, вероятно, является частью точки упражнения. Чтобы получить точное приближение, вы должны ожидать, что вам нужно, чтобы последний термин был небольшим, поэтому вам нужно n! для превышения x^n. По приближению Стирлинга n! ~ (n/e)^n, поэтому вы хотите n> e * x, где e = 2.71828 ..., поэтому n> = 33. В этот момент ошибка равна 0,005, а увеличение n на c уменьшает размер ошибки примерно на коэффициент e^c.
• Вы должны ожидать больших ошибок в арифметике с двойной точностью, когда вычитаете большие числа на пути к получению небольшого конечного результата. Это, вероятно, здесь не проблема, так как наибольшие члены имеют величину около 2^14. Вы по-прежнему получаете достаточную точность, чтобы не заметить, что вы не можете приблизиться к греху (12), добавив больше терминов.

Answer 2

Ваш факторный() возвращает в длинный тогда это возвращаемый тип двойной.

И самое важное, long не может содержать 16!. Какой двойной может.

Answer 3

Фактор между двумя членами ряда MacLauring или Тейлор для синуса

-x*x/(2*k*(2k+1)) 

Это может быть выгодно использован, чтобы избежать всех полномочий и факториалы и их переполнение.

mxx = -x*x; 
term = mxx/6; 
sum = 1+term; 
k=2; 
while(not good enough) 
    term = term*mxx/(2*k*(2*k+1)); 
    sum = sum + term; 
    k = k+1; 

return sum*x 

Answer 4

Есть несколько вещей, которые я бы сделал по-другому здесь.

Я немного предвзятый от необходимости писать численное программное обеспечение, которое не только дает правильные результаты, но и быстро их получает, но я вижу много расходомерных вычислений здесь. Рассмотрим два последовательных ненулевых члена ряда, например (x^13)/(13!) И - (x^15)/(15!). Если вы уже знаете значение (x^13)/(13!), Для чего вам нужно больше вычислить, чтобы получить - (x^15)/(15!)? Ответ, намного меньше, чем нужно, чтобы вычислить (x^13)/(13!) в первую очередь. Если вы рабски следуете обычной формуле для серии MacLaurin и пересчитайте факториал для каждого нового термина, , чтобы получить 15! вы повторите все вычисления, которые вы уже сделали за 13 !, а затем выполните только два новых умножения. Чтобы вычислить эту серию эффективно без потерь в расчете и без , вводя ненужные большие числа (и все возможные проблемы они могут вызвать, даже если вы используете с плавающей запятой для всего), просто посмотрите на соотношение между одним ненулевым членом и следующий. Это легко вычислить и не требует перетаскивания в функции pow. Хороший алгоритм позволит вам легко увеличить количество терминов до тех пор, пока последний не будет иметь порядок точности чисел с плавающей запятой, при условии, что вы используете разумное входное значение x. (Нет необходимости в том, чтобы x был больше 2 * pi, или, если на то пошло, действительно нет веской причины устанавливать x больше, чем pi/2, так как синус любого большего размера значение x можно найти по значениям . синусоидальной для ввода в этом диапазоне)

Answer 5

Вот как я это сделать в Java (легко преобразовать в C++):

package math; 

/** 
* Sine using Maclauren series 
* @author Michael 
* @link http://stackoverflow.com/questions/29445615/approximating-of-sine-and-the-remainder 
* @since 4/4/2015 8:37 PM 
*/ 
public class Sine { 

    private static final double TWO_PI = 2.0*Math.PI; 
    private static final int numPoints = 21; 
    private static final int numTerms = 21; 

    public static void main(String[] args) { 
     double x = -Math.PI; 
     double dx = 2.0*Math.PI/(numPoints-1); 
     for (int i = 0; i < numPoints; ++i) { 
      System.out.println(String.format("# terms: %d angle (radians) %10.6f sine: %15.10f", numTerms, x, sine(x, numTerms))); 
      x += dx; 
     } 
    } 

    public static double sine(double radians, int numTerms) { 
     double value = 0.0; 
     // Start by making sure the angle -pi/2 <= x <= +pi/2 
     double x = normalizeAngle(radians); 
     double term = x; 
     for (int n = 3; n < numTerms; n += 2) { 
      value += term; 
      term *= -x*x/n/(n-1); 
     } 
     return value; 
    } 

    public static double normalizeAngle(double radians) { 
     return radians - TWO_PI*Math.floor((radians+Math.PI)/TWO_PI); 
    } 
} 

Вот результат:

java math.Sine 
# terms: 21 angle (radians) -3.141593 sine: -0.0000000224 
# terms: 21 angle (radians) -2.827433 sine: -0.3090169974 
# terms: 21 angle (radians) -2.513274 sine: -0.5877852526 
# terms: 21 angle (radians) -2.199115 sine: -0.8090169944 
# terms: 21 angle (radians) -1.884956 sine: -0.9510565163 
# terms: 21 angle (radians) -1.570796 sine: -1.0000000000 
# terms: 21 angle (radians) -1.256637 sine: -0.9510565163 
# terms: 21 angle (radians) -0.942478 sine: -0.8090169944 
# terms: 21 angle (radians) -0.628319 sine: -0.5877852523 
# terms: 21 angle (radians) -0.314159 sine: -0.3090169944 
# terms: 21 angle (radians) 0.000000 sine: 0.0000000000 
# terms: 21 angle (radians) 0.314159 sine: 0.3090169944 
# terms: 21 angle (radians) 0.628319 sine: 0.5877852523 
# terms: 21 angle (radians) 0.942478 sine: 0.8090169944 
# terms: 21 angle (radians) 1.256637 sine: 0.9510565163 
# terms: 21 angle (radians) 1.570796 sine: 1.0000000000 
# terms: 21 angle (radians) 1.884956 sine: 0.9510565163 
# terms: 21 angle (radians) 2.199115 sine: 0.8090169944 
# terms: 21 angle (radians) 2.513274 sine: 0.5877852526 
# terms: 21 angle (radians) 2.827433 sine: 0.3090169974 
# terms: 21 angle (radians) 3.141593 sine: -0.0000000224 

Process finished with exit code 0