Физическая модель движения

 

 

1. Введение

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

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

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

Исходя из поставленной цели, можно выделить следующие задачи: рассмотреть физическую и математическую модели решения предложенного задания; разработать алгоритм решения задачи реализованный в виде блок-схемы; написать работоспособную программу на языке программирования Pascal.

 

 

2. Постановка задачи и определение конечных целей.

Задано уравнение движения материальной точки . Необходимо построить физическую и математическую модель движения.

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

Поэтому необходимо решить задачу Коши для дифференциального уравнения 1 – го порядка на отрезке [tо=1, tк=2] с шагом h=0.1  и  начальным условием x(tо)= 16.

Ответ должен  быть получен в виде таблицы результатов:

 

t	x(1)	x(2)
t0	x0(1)	x0(2)
t1	x1(1)	x1(2)
…	…	…
tk	xk(1)	xk(2)

 

где:  x(1), x(2)  - решения,  полученные различными  численными  методами.

Данные таблицы визуализировать в виде графиков.

Чтобы решить обыкновенное дифференциальное уравнение, необходимо знать значения зависимой переменной и (или) её производных при некоторых значениях независимой переменной. Если эти дополнительные условия задаются при одном значении независимой переменной, то такая задача называется задачей с начальными условиями, или задачей Коши.

Задачу Коши можно сформулировать следующим образом.

Пусть дано дифференциальное уравнение dx/dt=f(t,x)   и начальное условие x(t0) = x0.  Требуется найти функцию x(t), удовлетворяющую как указанному уравнению, так и начальному условию.

Численное решение задачи Коши сводится к табулированию искомой функции.

График решения дифференциального уравнения называется интегральной кривой.

Геометрический смысл задачи:

x’=f(t,x)  - тангенс угла наклона касательной к графику решения в точке (t, x) к оси ОХ, - угловой коэффициент (рис.1).

 

 

Существование решения:

Если правая часть f(t,x) непрерывна в некоторой области R, определяемой неравенствами

|t – t0| < a ;  |x – x0| < b

то существует, по меньшей мере, одно решение  x = x(t), определённое в окрестности  |t – x0| < h, где h – положительное число.

Это решение единственно, если в R выполнено условие Липшица

                             

где N – некоторая постоянная (константа Липшица), зависящая, в общем случае, от a и b. Если f(x, y)имеет ограниченную производную в R, то можно положить

N = max|

|    при   .

 

3. Разработка физической и математической модели задачи.

При использовании численных методов выполняется замена отрезка [t0, T] – области непрерывного изменения аргумента t множеством , состоящего из конечного числа точек  t0 < t1 < … < tn = T – сеткой.

При этом ti называют узлами сетки.

Во многих методах используются равномерные сетки с шагом

                                           

Задача Коши, определённая ранее на непрерывном отрезке [t0, T], заменяется её дискретным аналогом – системой уравнений, решая которую можно последовательно найти значения x1, x2, …, xn – приближённые значения функции в узлах сетки.

                                  (1)

Численное решение задачи Коши широко применяется в различных областях науки и техники, и число разработанных для него методов достаточно велико.

В данном случае будут использоваться одношаговые методы, в которых для нахождения следующей точки на кривой x = f(t) требуется информация лишь об одном предыдущем шаге. Одношаговыми являются метод Эйлера и методы Рунге – Кутта.

Применив вышеописанные знания на нашем примере, мы имеем:

дифференциальное уравнение первого порядка

x’ = f(t, x)

с начальным условием

x0(t0)= x0(1) = 16.

Выберем шаг h и введём обозначения:

                                   ti = t0 + i.h   и  xi = x(ti) ,

где    i = 0, 1, 2, …

         ti – узлы сетки,

        xi- значение интегральной функции в узлах .

 

Для решения поставленной задачи методом Эйлера модифицированным:

При использовании модифицированного метода Эйлера шаг h делится на два отрезка.

Иллюстрации к решению приведены на рисунке 1.

Рисунок 1 – Метод Эйлера модифицированный

 

 

Проведем решение в несколько этапов.

1. Обозначим точки: A(ti, xi), C(ti+h/2, xi+h/2*f(ti,xi)) и B(ti+1, xi+1).
2. Через точку  А проведем прямую под углом α, где

tgα = f(ti,xi)

3. На этой прямой найдем точку C(ti+h/2, xi+h/2*f(ti,xi)).
4. Через точку  С проведем прямую под углом α1, где

tgα1 = f(ti + h/2, xi + h/2 * f(ti,xi))

5. Через точку А проведем прямую, параллельную последней прямой.
6. Найдем точку В(ti+1, xi+1). Будем считать В(ti+1, xi+1) решением дифференциального уравнения при t=ti+1.
7. После проведения вычислений, получим формулу для определения значения xi+1:

xi+1 = xi + h* f(ti + h/2, xi + h/2 * f(ti,xi))                      

Модифицированный метод Эйлера дает меньшую погрешность. На рисунке 4 это хорошо видно. Так величина ε1 характеризует погрешность метода Эйлера, а ε – погрешность метода Эйлера модифицированного.

Блок-схема процедуры решения дифференциального уравнения методом Эйлера модифицированным приведена на рисунке 2.

 

F(t, x) – заданная функция – должна быть описана отдельно.

Входные параметры:

T0, TK – начальное и конечное значения независимой переменной;

X0 – значение x0 из начального условия  x(t0) = x0;

h – шаг;

Выходные параметры:

X – массив значений искомого решения в узлах сетки.

 

Для решения поставленной задачи методом Рунге – Кутта:

Согласно методу Рунге – Кутта четвёртого порядка, последовательные значения xi искомой функции x определяются по формуле:

xi+1 = xi + ∆xi                                            

где

∆xi = 1/6 * (k1 + 2*k2 + 2*k3 + k4),

i = 0, 1, 2, …   

а числа  k1(i),  k2(i),  k3(i),  k4(i)  на каждом шаге вычисляются по формулам:

 

k1 = h*f(ti, xi)

k2 = h*f(ti + h/2, xi + k1/2)

k3 = h*f(ti + h/2, xi + k2/2)

k4 = h*f(ti + h, xi + k3)

Это явный четырёхэтапный метод четвёртого порядка точности.

Методы Рунге – Кутта легко программируются и обладают значительной точностью и устойчивостью для широкого круга задач.

На рисунке 3 приведена блок-схема процедуры Runge(T0, TK, X0, h, X) для решения задачи Коши описанным выше методом Рунге – Кутта.

 

F(x, y) – заданная функция – должна быть описана отдельно.

Входные параметры:

T0, TK – начальное и конечное значения независимой переменной;

X0 – значение x0 из начального условия  x(t0) = x0;

h – шаг;

N – количество итераций.

Выходные параметры:

X – массив значений искомого решения в узлах сетки;

 

4. Разработка алгоритма в виде блок-схемы.

На рисунке 4 приведена блок-схема алгоритма основной программы для решения задачи Коши и получения результатов с фиксированным количеством отрезков разбиения N. В основной программе происходит обращение к процедурам Eiler(T0, TK, X0, h, X) и Runge (T0, TK, X0, h, X), вычисляющей значения искомой функции xj в точках tj методами Эйлера модифицированным и Рунге – Кутта.

Исходными данными в данной задаче являются:

T0, TK – начальное и конечное значения независимой переменной;

X0 – значение x0 из начального условия  x(t0) = x0;

h – шаг.

Результаты работы программы выводятся в виде трех столбцов:

T – массив значений узлов сетки;

X1 – массив значений искомого решения в соответствующих узлах сетки, полученный при решении методом Эйлера модифицированным;

X2 – массив значений искомого решения в соответствующих узлах сетки, полученный при решении методом Рунге-Кутта.

 

5. Разработка программы для ЭВМ.

Листинг программы:

 

program zadacha;

uses crt, graphABC;

var

    x1 : array [0..11] of real;

    x2 : array [0..11] of real;

    t  : array [0..11] of real;

    r, xnull, ynull, x, y, i, t0, tk, x0, N, a, b : integer;

    h: real;

 

procedure Eiler(t0,tk,x0,N:integer; h:real; var x : array [0..11] of real);

var

    f1, f2 : real;

begin

    t[0]  := t0;

    t[N]  := tk;

    x[0] := x0;

 

    for i:= 0 to N do

    begin

        t[i] := t[0] + i*h;

        f1  := 2*ln(t[i])*sqrt(x[i]);

        f2   := 2*ln(t[i]+h/2)*sqrt(x[i]+h/2*f1);

 

        x[i+1] := x[i] + h*f2;

 

    end;

end;

 

procedure Runge(t0,tk,x0,N:integer; h:real; var x : array [0..11] of real);

var

    k1, k2, k3, k4, f, k : real;

begin

    t[0]  := t0;

    t[N]  := tk;

    x[0] := x0;

 

    for i:= 0 to N do

    begin

        t[i] := t[0] + i*h;

        f    := 2*ln(t[i])*sqrt(x[i]);

        k1   := h*f;

        f    := 2*ln(t[i]+h/2)*sqrt(x[i] + k1/2);

        k2   := h*f;

        f    := 2*ln(t[i]+h/2)*sqrt(x[i] + k2/2);

        k3   := h*f;

        f    := 2*ln(t[i]+h)*sqrt(x[i] + k3);

        k4   := h*f;

        k    := (k1 + 2*k2 + 2*k3 + k4)/6;

 

        x[i+1] := x[i]+k;

    end;

end;

 

begin

    clrscr;

    // входные параметры //

    t0 := 1;

    tk := 2;

    x0 := 16;

    h  := 0.1;

    N  := 10;

    // функция x:= 2*ln(t)*sqrt(x) //

  

    Eiler (t0,tk,x0,N,h,x1);

    Runge (t0,tk,x0,N,h,x2);

 

    writeln(‘Таблица результатов:’);

    write(‘t’:9);

    write(‘x1’:9);

    writeln(‘x2’:10);

    for i:= 0 to N do

    begin

        write(t[i]:10:2);

        write(x1[i]:10:2);

        writeln(x2[i]:10:2);

    end;

    readln();

 

    xnull := 290;  // начало координат

    ynull := 220;

    Line(xnull,ynull,xnull+300,ynull);

    Line(xnull,ynull-214,xnull,ynull);

    Line(xnull,ynull,xnull-30,ynull);

    Line(xnull,ynull,xnull,ynull+30);

    Line(xnull+300,ynull,xnull+290,ynull-7);

    Line(xnull+300,ynull,xnull+290,ynull+7);

    Line(xnull,ynull-214,xnull+7,ynull-204);

    Line(xnull,ynull-214,xnull-7,ynull-204);

    TextOut(xnull+300,ynull+10, ‘t’);

    TextOut(xnull-23,ynull-214, ‘x’);

 

    r := 3;

    SetPenColor(clRed);

    for i:= 0 to (N-1) do

    begin

        x := round(t[i]*10)*10;

        y := round(x1[i]*10);

        MoveTo(x+xnull,ynull-y);

        circle(x+xnull,ynull-y,r);

        a := round(t[i+1]*10)*10;

        b := round(x1[i+1]*10);

        LineTo(a+xnull, ynull-b);

        circle(a+xnull,ynull-b,r);

    end;

 

    Line(xnull+240,ynull-50,xnull+270,ynull-50);

    TextOut(xnull+280,ynull-60, ‘x1’);

 

 

        r := 3;

 

    for i:= 0 to (N-1) do

    begin

        SetPenColor(clRed);

        x := round(t[i]*10)*10;

        y := round(x1[i]*10);

        MoveTo(x+xnull,ynull-y);

        circle(x+xnull,ynull-y,r);

        a := round(t[i+1]*10)*10;

        b := round(x1[i+1]*10);

        LineTo(a+xnull, ynull-b);

        circle(a+xnull,ynull-b,r);

        if (i=0) then

        begin

              SetPenColor(clBlack);

              circle(x+xnull,ynull, r);

              circle(xnull,ynull-y, r);

              TextOut(x+xnull-7,ynull+5, ‘t0’);

              TextOut(xnull-21,ynull-y-7, ‘x0’);

        end;

        if (i=(N-1)) then

        begin

              SetPenColor(clBlack);

              circle(x+xnull,ynull, r);

              circle(xnull,ynull-y, r);

              TextOut(x+xnull-7,ynull+5, ‘tk’);

              TextOut(xnull-21,ynull-y-7, ‘xk’);

        end;

    end;

    SetPenColor(clRed);

    Line(xnull+240,ynull-50,xnull+270,ynull-50);

    TextOut(xnull+280,ynull-60, ‘x1’);

     readln();

 

    r := 3;

    for i:= 0 to (N-1) do

    begin

        SetPenColor(clGreen);

        x := round(t[i]*10)*10;

        y := round(x2[i]*10);

        MoveTo(x+xnull,ynull-y);

        circle(x+xnull,ynull-y,r);

        a := round(t[i+1]*10)*10;

       b := round(x2[i+1]*10);

        LineTo(a+xnull, ynull-b);

        circle(a+xnull,ynull-b,r);

        if (i=0) then

        begin

              SetPenColor(clBlack);

              circle(x+xnull,ynull, r);

              circle(xnull,ynull-y, r);

              TextOut(x+xnull-7,ynull+5, ‘t0’);

              TextOut(xnull-21,ynull-y-7, ‘x0’);

        end;

        if (i=(N-1)) then

        begin

              SetPenColor(clBlack);

              circle(x+xnull,ynull, r);

              circle(xnull,ynull-y, r);

              TextOut(x+xnull-7,ynull+5, ‘tk’);

              TextOut(xnull-21,ynull-y-7, ‘xk’);

        end;

    end;

 

    SetPenColor(clGreen);

    Line(xnull+240,ynull-30,xnull+270,ynull-30);

    TextOut(xnull+280,ynull-40, ‘x2’);

 

end.

 

6. Работа по программе и обработка результатов.