Код Matlab для аппроксимации экспоненциальной функции

Question

Кто-нибудь знает, как сделать так, чтобы следующий код Matlab более точно аппроксимировал экспоненциальную функцию при работе с большими и отрицательными действительными числами?

Например, когда x = 1, код работает хорошо, когда x = -100, он возвращает ответ 8.7364e+31, когда он должен быть ближе к 3.7201e-44.

Код выглядит следующим образом:

s=1
a=1;
y=1;
for k=1:40
    a=a/k;
    y=y*x;
    s=s+a*y;
end
s

Любая помощь приветствуется, ура.

РЕДАКТИРОВАТЬ:

Итак, вопрос в следующем:

Какую математическую функцию аппроксимирует этот код? (Я говорю экспоненциальная функция.) Работает ли она, когда x = 1? (Да.) К сожалению, использование этого, когда x = -100 дает ответ s = 8.7364e+31. Ваш коллега считает, что в программе есть глупая ошибка, и просит вашей помощи. Тщательно объясните поведение и дайте простое исправление, которое даст лучший результат. [Вы должны предложить модификацию вышеприведенного кода или его использование. Вы также должны проверить свои простые исправления.]

Поэтому я в некоторой степени понимаю, что проблема связана с большими числами, когда между терминами 16 (или более) порядков, точность теряется, но решение ускользает от меня.

Спасибо

РЕДАКТИРОВАТЬ:

Итак, в конце концов я пошел с этим:

s = 1;
x = -100;
a = 1;
y = 1;
x1 = 1;

for k=1:40
    x1 = x/10;
    a = a/k;
    y = y*x1;
    s = s + a*y;
end

s = s^10;
s

Не уверен, что это полностью правильно, но возвращает хорошие приближения.

exp (-100) = 3,720075976020836e-044
s = 3,722053303838800e-044

После дальнейшего анализа (и, к сожалению, отправки задания) я понял, что увеличение количества итераций и, следовательно, увеличение сроков еще более повышает эффективность. На самом деле следующее было еще более эффективным:

s = 1;
x = -100;
a = 1;
y = 1;
x1 = 1;

for k=1:200
    x1 = x/200;
    a = a/k;
    y = y*x1;
    s = s + a*y;
end

s = s^200;
s

Который дает:

exp (-100) = 3,720075976020836e-044
s = 3,720075976020701e-044

Answer 1

Как отмечает Джон в комментарии, внутри цикла есть ошибка. Строка y = y*k не делает то, что вам нужно. Посмотрите внимательнее на условия в серии для exp(x).

В любом случае, я предполагаю, что именно поэтому вам было дано задание по домашнему заданию, чтобы узнать, что подобные серии не очень хорошо сходятся для больших значений. Вместо этого вы должны подумать, как сделать уменьшение диапазона.

Например, вы можете использовать личность

exp(x+y) = exp(x)*exp(y)

в ваших интересах? Предположим, вы храните значение exp(1) = 2.7182818284590452353...

Теперь, если бы я попросил вас вычислить значение exp(1.3), как бы вы использовали вышеуказанную информацию?

exp(1.3) = exp(1)*exp(0.3)

Но мы уже ЗНАЕМ значение exp (1). На самом деле, если немного подумать, это позволит вам уменьшить диапазон экспоненты до того, чтобы ряды быстро сходились только для abs(x) <= 0.5.

Редактировать: есть второй способ уменьшить диапазон, используя вариацию одного и того же персонажа.

exp(x) = exp(x/2)*exp(x/2) = exp(x/2)^2

Таким образом, предположим, что вы хотите вычислить экспоненту большого числа, возможно, 12,8. Чтобы добиться этого достаточно быстро, чтобы сойтись, потребуется много терминов в простых сериях, и произойдет много субтрактивного отмены, так что вы все равно не получите хорошую точность. Тем не менее, если мы признаем, что

12.8 = 2*6.4 = 2*2*3.2 = ... = 16*0.8

тогда, если бы вы могли эффективно вычислить экспоненту 0,8, тогда желаемое значение легко восстановить, возможно, путем повторного возведения в квадрат.

exp(12.8)
ans =
          362217.449611248

a = exp(0.8)
a =
          2.22554092849247
a = a*a;
a = a*a;
a = a*a;
a = a*a
          362217.449611249

exp(0.8)^16
ans =
          362217.449611249

Обратите внимание, что всякий раз, когда вы выполняете сокращение диапазона с помощью подобных методов, в то время как вы можете столкнуться с числовыми проблемами из-за необходимых дополнительных вычислений, вы обычно выходите далеко вперед из-за значительно улучшенной конвергенции ваших рядов.

Answer 2

Почему вы думаете, что это неправильный ответ? Посмотрите на последний член этой последовательности и его размер, и скажите мне, почему вы ожидаете, что у вас должен быть ответ, близкий к 0.

В моем первоначальном ответе говорилось, что проблема заключалась в ошибке округления. Это будет проблемой для этого базового подхода, но почему вы думаете, что 40 достаточно терминов для соответствующего математического (в отличие от компьютерной арифметики с плавающей запятой) ответа.

100 ^ 40/40! ~ = 10 ^ 31.

Woodchip имеет правильную идею с уменьшением диапазона. Это типичный подход, который люди используют для очень быстрой реализации этих функций. Как только вы все это выяснили, вы имеете дело с ошибками округления чередующихся последовательностей путем суммирования смежных членов в цикле и пошагового выполнения с k = 1: 2: 40 (например). Это не сработает, пока вы не воспользуетесь идеей крошки, потому что для x = -100 слагаемые растут очень долго. Вам нужно | х | < 1, чтобы гарантировать, что промежуточные сроки сокращаются, и, таким образом, перезапись будет работать.