Метод Эйлера - Весенние колебания в MatLab

Question

Если у меня есть дифференциальное уравнение второго порядка в терминах дифференциала первого порядка, описывающего систему «масса-пружина-демпфер», как я могу использовать метод Эйлера для построения этого уравнения, когда я не знаю первое-дифференциал?Я хочу сделать это в MatLab.Это домашнее задание, поэтому я не написал ни одного кода ... Я просто хочу кратко описать, как вы должны это делать.

Производная второго порядка: d2 (t + 1) = (-1 / m) * (c * d1 + k * y), где c, m, k - постоянные, y изначально равно 1, а d1это дифференциал первого порядка, который начинается с 0, а t - это время.

Есть идеи?

Спасибо :)).

Answer 1

Уравнение 2-го порядка можно преобразовать в систему дифференциальных уравнений первого порядка.

function dy = ex(y)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -c/m*y(2) - k/m*y(1);

Отсюда вы можете использовать встроенный решатель Matlab.ode23s будет работать просто отлично:

[t,y] = ode23s(@ex, y0, tspan)

Answer 2

Перепишите уравнение второго порядка как систему уравнений первого порядка. Неизвестные соответствуют положению и скорости.

Answer 3

Система дифференциальных уравнений первого порядка может выглядеть следующим образом, тогда как y1 = y. ' используется для обозначения производной по времени.

y1' = y2
y2' = -c/m*y2 - k/m*y1

Answer 4

Существует аналитическое решение вида (x (t), x '(t)) = exp (A t) * (x (0), x' (0)), где A - матрица 2x2.Если вам не нужно использовать ODE-решатель Matlab, это верный способ найти временную эволюцию системы.

Чтобы найти A, напишите свою систему в этой форме

• X = (x, x ')
• dX = A * X
• X = expm (A t) X (0)% обратите внимание, что для взятия экспоненты матрицы необходимоиспользовать expm, иначе вы просто получите экспоненту каждого элемента в отдельности

, так что здесь A = [0 1;-k / m -c / m]

установка k / m = 1 и c / m = 0,1, мы можем написать

t = linspace(0,20, 1000);
A = [0 1; -1 -0.1];
for j = 1:length(t); x(j,:) = expm([0 1; -1 -0.1]*t(j))*[1;0]; end
plot (t,x)
legend ('position', 'velocity');
title ('underdamped spring starting at y = 1; y'' = 0')