Составление  уравнений  движения

Уравнения высших порядков приводят к системе уравнений первого порядка введением дополнительных переменных. Например, для уравнения второго порядка

y" = f(x,y',y)

примем v= y'. Тогда y" = v' и имеем систему уравнений:

v’ = f(x, y, v);

y’ = v.

Графическая интерпретация численного решения обыкновенного дифференциального уравнения показана на рис. 7 на примере простейшего метода Эйлера. Известной является функция y0 в точке x0.

Рис. 7. Графическая интерпретация метода Эйлера

Решение находится для ряда значений независимой переменной х с шагом h:

x1 = x0 + h ; x2 = x1 + h ; . xn+1 = xn + h .

Значение y1 (рис. 7) находится на пересечении прямой, проведенной из точки (х0,у0) под углом a0 = arctg(y0') и перпендикуляра, проведенного к оси абсцисс из точки х1. Процесс последовательно повторяется для других зна­чений х:

y1 = y0 + h×y0' = y0 + h×f(x0,y0) ,

y2 = y1 + h×y1' = y1 + h×f(x1,y1), .

yn+1 = yn + h×yn' = yn + h×f(xn,yn) .

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

yn+1 = yn + h·yn' +h2/2·y” +h3/6·y(3) + .

Метод Эйлера учитывает производную только первого порядка, поэтому является методом первого порядка. Чаще всего используется метод Рунге-Кутта четвертого поряд­ка, алгоритм которого имеет вид:

yn+1 = yn + (k1 + 2×k2 +2×k3 + k4) / 6 ,

где k1 = h×f(xn,yn) ;

k2 = h×f(xn + 0.5×h,yn + 0.5×k1) ;

k3 = h×f(xn + 0.5×h,yn + 0.5×k2) ;

k4 = h×f(xn + h,yn + k3) .

Начало программы (вариант) для решения дифференциального уравнения второго по­рядка

m×x" + b×x' + c×x = F

показано ниже.

Program DIFUR;

Uses Crt,Dos,Lib,Graph;

type mas = array[1 10] of real;

Var Rez :text; d,y,y1 :mas;

Filerez,text,Sxmax,Sxst,Stmax :string[50];

m,b,c,F,h,hp,t,tp,tmax,xmax,xst,nx,ny,nmax :real;

n :integer;

KL,nd,j :byte;

Procedure Prav(var y,d:mas);

í y[1] = x' - скорость массы d[1] = x''- ускорение массы ý

í y[2] = x - перемещение массы d[2] = x' - скорость массы ý

begin

d[1] := (F - b*y[1] - c*y[2])/m; d[2] := y[1];

end;

Procedure Rynge;

Var j :byte; yy,k :array[1 nd] of real;

begin

Prav(y,d);

For j:=1 to nd do begin

yy[j]:=y[j]; k[j]:=h*d[j]; y1[j]:=yy[j] + 0.5*k[j]; y[j]:=y[j] + k[j]/6; end;

t:=t+0.5*h; Prav(y1,d);

For j:=1 to nd do begin

k[j]:=h*d[j]; y1[j]:=yy[j] + 0.5*k[j]; y[j]:=y[j] + k[j]/3; end;

Prav(y1,d);

For j:=1 to nd do begin

k[j]:=h*d[j]; y1[j]:=yy[j] + k[j]; y[j]:=y[j] + k[j]/3; end;

t:=t+0.5*h; Prav(y1,d);

For j:=1 to nd do begin k[j]:=h*d[j]; y[j]:=y[j] + k[j]/6; end;

end;

Procedure Start; .

Процедура Prav (var y,d:max) предназначена для вычисления пра­вых частей уравнений. Использованные в ней массивы типа mas: y – вектор - решение (выходные параметры); d - производные век­тор - решения (производные выходных параметров).

Процедура Rynge реализует метод Рунге-Кутта четвертого поряд­ка.