Matlab Дифференциальные уравнения Метод Эйлера

Мне нужна помощь в составлении дифференциального уравнения... оно все время выходит странным, и график не такой, каким он должен быть.

function [dydt] = diff(y,t)

dydt = (-3*y)+(t*(exp(-3*t)));

end

tI = 0;
yI = -0.1;
tEnd = 5;
dt = 0.5;

t = tI:dt:tEnd;
y = zeros(size(t));
y(1) = yI;

for k = 2:numel(y)
    yPrime = diff(t(k-1),y(k-1));
    y(k) = y(k-1) + dt*yPrime;
end

plot(t,y)
grid on
title('Engr')
xlabel('Time')
ylabel('y(t)')
legend(['dt = ' num2str(dt)])

Это мой код, но график совсем не похож на то, на что он должен быть похож. Я пропустил что-то вроде индекса для оператора for?

Изменить

Я получаю сообщение об ошибке:

Error using diff
Difference order N must be a positive integer scalar.

Error in diff3 (line 12)
    yPrime = diff(t(k-1),y(k-1));

matlab differential-equations
person Apetro    schedule 06.07.2013    source источник
comment
У вас проблемы с именами функций: в Matlab уже есть встроенная функция 'diff' и попробуйте вызвать эту функцию вместо вашей собственной функции. Попробуйте использовать другое имя для вашей функции «diff».   -  person Danil Asotsky    schedule 06.07.2013
comment
Другая ваша проблема заключается в том, что ваша функция ODE определена как diff(y,t), но вы звоните с перевернутыми аргументами: diff(t(k-1),y(k-1)). Я бы рекомендовал следовать тому, что делают решатели ОДУ Matlab (например, ode45), и определять аргументы вашей функции ОДУ как t, а затем y. Вы также обнаружите, что ваш временной шаг dt слишком велик даже для этой простой функции. Попробуйте 0,1 или меньше.   -  person horchler    schedule 06.07.2013
comment
Вы также захотите, чтобы ваша функция ODE была подфункцией (или отдельной функцией M-файла), а не началом вашего кода, как показано здесь.   -  person horchler    schedule 06.07.2013

Ответы (1)


arrow_upward
1
arrow_downward

После исправления ошибок, на которые указали Данил Асоцкий и horchler в комментариях:

  1. избежание конфликта имен со встроенной функцией 'diff'
  2. изменение порядка аргументов на t,y.
  3. уменьшение временного шага dt до 0,1
  4. преобразование правой части ОДУ в анонимную функцию

(и удалив ненужные скобки в определении функции), ваш код может выглядеть так:

F = @(t,y) -3*y+t*exp(-3*t);

tI = 0;
yI = -0.1;
tEnd = 5;
dt = 0.1;

t = tI:dt:tEnd;
y = zeros(size(t));
y(1) = yI;

for k = 2:numel(y)
    yPrime = F(t(k-1),y(k-1));
    y(k) = y(k-1) + dt*yPrime;
end

plot(t,y)
grid on
title('Engr')
xlabel('Time')
ylabel('y(t)')
legend(['dt = ' num2str(dt)])

который работает так, как ожидалось:

output

person Community    schedule 05.01.2015