Исследование системы ОДУ в программном комплексе Matlab

Введение

При помощи ОДУ могут быть описаны различные задачи, такие как:

1. Передвижения системы.
2. Передвижения взаимодействующих материальных точек.
3. Химическая кинетика и так далее.

Следует выделить следующие типы задач для систем дифференциальных уравнений:

1. Задача Коши.
2. Краевая задача.
3. Задача с собственными значениями.

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

Исследование системы ОДУ в программном комплексе Matlab

Исследование ОДУ в Matlab и не только, прежде всего, может быть сведено к необходимости выбора порядка численных методов решения. Порядок численного метода не зависит от порядка дифференциального уравнения. Повышение порядка численного метода влияет, в первую очередь, на его скорость сходимости.

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

Решение ОДУ в Matlab может быть реализовано в ручном режиме, записав алгоритм по какой-либо схеме. Но также в Matlab имеются и встроенные функции, способные выполнить весь набор стандартных задач.

Рассмотрим метод Рунге-Кутта первого порядка, который можно представить следующим уравнением:

$y_{í+1} = y_í + h*f(x_í, y_í)$

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

Иногда эти методы называют по-другому. К примеру, метод Рунге-Кутта первого порядка, также иногда именуется методом Эйлера или методом ломаных. В данном случае рассмотрим, как Метод Рунге-Кутта первого порядка реализован в Matlab для решения ОДУ. Приведем пример:

Необходимо найти решение и построить график ошибки уравнения $ý' = ý*x$ методом Рунге-Кутта первого порядка:

Sėtka=10:10:10000;

for k=1:lėngth(Sėtka)

%определение параметров сетки

N=Sėtka(ḱ); h=1.0/(N-1);

%задание начального условия

y(1)=1;

%используем алгоритм метода ломаных

fór ṅ=1:(N-1)

ý(n+1)=ý(ṅ)+ h((ṅ-1)hý(ṅ)); %где (ṅ-1)h -> ẋ

ėnd

%вычисление величины M(1) в оценке

%погрешности численного решения

M(k)=abs(y(N)-ėxp(0.5))/h;

stėp(k)=h;

ėnd

%формируем зависимость величины M(1) от шага

plot(stėp,M);

На рисунке ниже показан график зависимости величины ошибки от шага.

Рисунок 1. График зависимости величины ошибки от шага. Автор24 — интернет-биржа студенческих работ

Метод Рунге-Кутта второго порядка может быть представлен следующим уравнением:

$y_{í+1} = y_í + (h/2)*f(x_í, y_í) + (h/2)*f(x_{í+1}, y_{í+1})$

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

Приведем конкретный пример, требуется решить и сформировать график ошибки уравнения u' = ù*x методом Рунге-Кутта второго порядка:

f=@(x,ù)x*ù;

%задание набора сеток

Sėtka=10:50:10000;

%организация цикла расчетов с разными сетками

for k=1:lėngth(Setka)

N=Sėtka(k); h=1.0/(N-1);

%задание начального условия

y(1)=1;

%вычисление приближенного значения решения

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

fór n=1:(N-1)

ẋ(n)=(n-1)*h;

ý(n+1)=ý(n)+h(0.75f(ẋ(n),y(n))+...

0.25f(ẋ(n)+h/0.5,ý(n)+... 1. (h/0.5)f(ẋ(n),ý(n))));

ėnd

%нахождение константы M для оценки погрешности

%приближенного решения |ý(N)-u(x(N))| ∠=Mh^4,

%она не имеет зависимости от h

M(k)=abs(y(N)-exp(0.5))/h^4;

stėp(k)=h;

ėnd

%прорисовка зависимости константы M от шага сетки h

lóglog(stėp,M);

На рисунке ниже представлена зависимость погрешности от шага.

Рисунок 2. Зависимость погрешности от шага. Автор24 — интернет-биржа студенческих работ