m y'' + k2 y' + k1 y = f(t)

[logarithm]

Есть задача.

Твердое тело подвешено на пружине с амортизатором. На это тело действует вертикальная сила f(t) и его движение в вертикальной плоскости описывается следующим уравнением
m y'' + k2 y' + k1 y = f(t) где y – вертикальная координата тела, y' – производная y по по времени, y'' — вторая производная y по времени, m — масса тела, k2 – коэффициент амортизации.
С помощью языка Си++ проимитируйте 60 с движения тела с целью получения графика движения тела по вертикали при m = 1; k=0,1; k2=0,3; f(t) =1,0; при любом t.

Нужен примерный алгоритм.

-- Вт дек 21, 2010 18:19:30 --

Хотя бы подскажите в какую сторону копать.

[gris]

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

[Portnov]

Мнээ... Так при постоянном f(t) это уравнение же аналитически решается :) Решите на бумажке, получите формулу для y(t) и считайте по ней :)

[Circiter]

2logarithm
Есть простой (но практически непригодный) способ численного интегрирования таких уравнений. Сначала нужно выразить $y''=(f(t)-k_2 y'-k_1 y)/m$. Понятно, что $y'$ теперь можно получить просто проинтегрировав $y''$, например, по-Эйлеру; аналогично получается искомая $y$.

Псевдокод:

dt=0.1; // "Infinitesimal" time interval.
t=0;
// Your initial conditions here.

while(t<=60)
{
    ddy=(f(t)-k2*dy-k1*y)/m; // Calculate second derivative of y.

    dy+=ddy*dt; // Integrate it using Euler method.
    y+=dy*dt; // In the same way, get y.
    // Plot(t, y);

    t+=dt; // Update the time position.
}

А вообще, гляньте в какую-нибудь книжку по численным методам, подыщите что-нибудь вроде метода Рунге-Кутта... Можно посмотреть для начала хотя-бы в H.Press et al., Numerical Recipes in C, там упрощенное изложение.