Помощь по коду MATLAB. Обратный метод Эйлера

Question

Вот код MATLAB / FreeMat , который я получил для численного решения ODE с использованием метода обратного Эйлера . Тем не менее, результаты несовместимы с результатами моего учебника, а иногда даже смешно противоречивы. Что не так с кодом?

function [x,y] = backEuler(f,xinit,yinit,xfinal,h)

    %f - this is your y prime
    %xinit - initial X
    %yinit - initial Y
    %xfinal - final X
    %h - step size

    n = (xfinal-xinit)/h; %Calculate steps

    %Inititialize arrays...
    %The first elements take xinit and yinit corespondigly, the rest fill with 0s.
    x = [xinit zeros(1,n)];
    y = [yinit zeros(1,n)];

    %Numeric routine
    for i = 1:n
        x(i+1) = x(i)+h;
        ynew = y(i)+h*(f(x(i),y(i)));
        y(i+1) = y(i)+h*f(x(i+1),ynew);
    end
end

Answer 1

Ваш метод - это метод нового типа . Это ни назад, ни вперед, Эйлер. : -)

Форвард Эйлер: y1 = y0 + h*f(x0,y0)

Обратный Эйлер solve in y1: y1 - h*f(x1,y1) = y0

Ваш метод: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Ваш метод не отсталый Эйлер.

Вы не решаете в y1, вы просто оцениваете y1 с помощью прямого метода Эйлера. Я не хочу продолжать анализ вашего метода, но я считаю, что он действительно будет вести себя плохо, даже по сравнению с форвардом Эйлером, поскольку вы оцениваете функцию f в неправильной точке.

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

y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))

Answer 2

Посмотрите на числовые рецепты , в частности главу 16, интегрирование обыкновенных дифференциальных уравнений. Известно, что метод Эйлера имеет проблемы:

Существует несколько причин, по которым метод Эйлера не рекомендуется для практического использования, среди них: (i) метод не очень точен по сравнению с другими, более привлекательными, методами, работающими с эквивалентным размером шага, и (ii) также не является очень стабильный

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

Answer 3

Единственная проблема, которую я могу заметить, это то, что строка:

n=(xfinal-xinit)/h

Должно быть:

n = abs((xfinal-xinit)/h)

Чтобы избежать отрицательного количества шагов. Если вы двигаетесь в направлении отрицательного-x, убедитесь, что для функции задан отрицательный размер шага.

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

PS. Это не «обратный метод Эйлера», это просто обычный старый метод Эйлера.

Если это домашнее задание, пожалуйста, пометьте его так.

Answer 4

Если вы действительно не хотите решить ODE с помощью метода Эйлера, который вы написали самостоятельно, вам следует взглянуть на встроенные решатели ODE .

На sidenote: вам не нужно создавать x(i) внутри цикла, как это: x(i+1) = x(i)+h;. Вместо этого вы можете просто написать x = xinit:h:xfinal;. Также вы можете написать n = round(xfinal-xinit)/h);, чтобы избежать предупреждений.

---

Вот решатели, реализованные в MATLAB.

ode45 основан на явном Рунге-Кутта (4,5) формула, Пара Дорман-Принц. Это один шаг решатель - в вычислениях y (tn) он нужен только решение сразу предыдущий момент времени, y (tn-1). В Вообще, ode45 - лучшая функция для применить в качестве первой попытки для большинства проблемы.

ode23 является реализацией явная пара Рунге-Кутта (2,3) Богацки и Шампин. Может быть больше эффективнее, чем ode45 на сырой допуски и при наличии умеренная жесткость. Как ode45, ode23 одноступенчатый решатель

ode113 - это переменный порядок Адамс-Башфорт-Моултон PECE решатель. Это может быть более эффективным, чем ode45 в строгие допуски и когда ODE функция файла особенно дорого оценивать. ode113 является многоступенчатый решатель - обычно он нужен решения на несколько предыдущих моменты времени для вычисления текущего решение.

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

ode15s - решатель переменного порядка на основе численного дифференцирования формулы (NDF). По желанию, он использует формулы обратной дифференциации (BDFs, также известный как метод Гира) которые обычно менее эффективны. подобно ode113, ode15s - это многошаговый решатель. Попробуйте ode15s, если ode45 не работает или очень неэффективно, и вы подозреваете, что проблема жесткая, или при решении дифференциально-алгебраическая задача.

ode23s основан на модифицированной Формула Розенброка порядка 2. Потому что это одностадийный решатель, это может быть более эффективный, чем ode15s на сыром допуски. Это может решить некоторые виды Сложные проблемы, для которых ode15s не эффективный.

ode23t является реализацией Трапециевидное правило с использованием «свободного» интерполянт. Используйте этот решатель, если проблема только умеренно жесткая и вам нужно решение без числовых демпфирования. ode23t может решить DAE.

ode23tb - это реализация TR-BDF2, неявный Рунге-Кутта формула с первой стадией, которая является шаг трапециевидного правила и второй этап, который является отсталым формула дифференцирования второго порядка. По построению та же итерация матрица используется при оценке как этапы. Как и ode23s, этот решатель может быть более эффективным, чем ode15s в сырой допуски.

Answer 5

Код в порядке.Просто вы должны добавить еще один цикл внутри цикла.Чтобы проверить уровень последовательности.

if abs((y(i+1) - ynew)/ynew) > 0.0000000001 
    ynew = y(i+1);
    y(i+1) = y(i)+h*f(x(i+1),ynew);
end

Я проверил на наличие фиктивной функции, и результаты были многообещающими.

Answer 6

Я думаю, что этот код может работать. Попробуйте это.

for i =1:n
    t(i +1)=t(i )+dt;
    y(i+1)=solve('y(i+1)=y(i)+dt*f(t(i+1),y(i+1)');
    end