Численные методы интегрирования дифференциальных уравнений первого порядка

Размещено на http://www.allbest.ru/

Министерство образования и науки Российской Федерации

Калужский филиал федерального государственного бюджетного образовательного учреждения высшего профессионального образования

«Московский государственный технический университет имени Н.Э. Баумана»

ЛАБОРАТОРНАЯ РАБОТА №1

ДИСЦИПЛИНА: "Вычислительные методы теории управления"

ТЕМА: "Численные методы интегрирования дифференциальных уравнений первого порядка"

Выполнил:

Аксенов Н.В.

Калуга, 2014 г.

Цель работы: Изучение численных методов интегрирования дифференциальных уравнений первого порядка. В ходе лабораторной работы выполняются исследования различных методов интегрирования дифференциальных уравнений по точности вычисления и быстродействию построенных на их основе алгоритмов.

Вариант 1:

Рассмотрим решение данного уравнения на промежутке [2,10] с нулевыми начальными условиями.

Решим данное уравнение с помощью встроенного решателя ode45:

fun_1(1-2*t*x)/(t*(t-1));

[t,y]=ode45(@fun_1,[2 10],[0]); plot(t,y)

Рисунок 1. График решения, полученного с помощью встроенного решателя.

дифференциальный уравнение алгоритм программа

Рисунок 2. График решений, полученных численными методами с шагом 0.1.

Рисунок 3. График решений, полученных численными методами с шагом 0.01.

1. Метод Эйлера.

Рисунок 4. Структурная схема алгоритма метода Эйлера.

Листинг программы MATLAB реализации метода Эйлера:

function [t,x] = eilera_difur(funkc,interval,inital,h)

funkc=inline(funkc);

i=1;

x(1)=inital;

t(1)=interval(1);

while t(i)<=interval(2)

x(i+1)=x(i)+h*funkc(t(i),x(i));

t(i+1)=t(i)+h;

i=i+1;

end

2. Метод Эйлера-Коши.

Рисунок 5. Структурная схема алгоритма метода Эйлера-Коши.

Листинг программы MATLAB реализации метода Эйлера-Коши:

function [t,x] = eilera_koshi_difur(funkc,interval,inital,h)

funkc=inline(funkc);

i=1;

x(1)=inital;

t(1)=interval(1);

while t(i)<=interval(2)

x(i+1)=x(i)+h*funkc(t(i),x(i));

t(i+1)=t(i)+h;

x(i+1)=x(i)+h*(funkc(t(i),x(i))+funkc(t(i+1),x(i+1)))/2;

i=i+1;

end

3. Метод Рунге-Кутта 3 и 4 порядков.

Рисунок 6. Структурные схемы алгоритмов метода Рунге-Кутта 3 и 4 порядков.

Листинг программы MATLAB реализации метода Рунге-Кутта 3 порядка:

function [t,x] = Rungekutta3(funkc,interval,inital,h)

funkc=inline(funkc);

i=1;

x(1)=inital;

t(1)=interval(1);

while t(i)<=interval(2)

K1=h*funkc(t(i),x(i));

K2=h*funkc(t(i)+h/2,x(i)+K1/2);

K3=h*funkc(t(i)+h,x(i)-K1+2*K2);

x(i+1)=x(i)+(K1+4*K2+K3)/6;

t(i+1)=t(i)+h;

i=i+1;

end

Листинг программы MATLAB реализации метода Рунге-Кутта 4 порядка:

function [t,x] = Rungekutta4(funkc,interval,inital,h)

funkc=inline(funkc);

i=1;

x(1)=inital;

t(1)=interval(1);

while t(i)<=interval(2)

K1=h*funkc(t(i),x(i));

K2=h*funkc(t(i)+h/2,x(i)+K1/2);

K3=h*funkc(t(i)+h/2,x(i)+K2/2);

K4=h*funkc(t(i)+h,x(i)+K3);

x(i+1)=x(i)+(K1+2*K2+2*K3+K4)/6;

t(i+1)=t(i)+h;

i=i+1;

end

4. Экстраполяционный метод Адамса 3 порядка.

Рисунок 7. Структурная схема алгоритма метода интерполяционного метода Адамса 4 порядка.

Листинг программы MATLAB реализации интерполяционного метода Адамса 4 порядка.

function [t,x] = Adamsa4(funkc,interval,inital,h)

t(1)=interval(1);

t(2)=t(1)+h;

t(3)=t(2)+h;

t(4)=t(3)+h;

[bt,bx]= eilera_koshi_difur(funkc,[interval(1) interval(1)+3*h],inital,h/10);

x(1)=inital;

x(2)=bx(11);

x(3)=bx(21);

x(4)=bx(31);

funkc=inline(funkc);

i=4;

while t(i)<interval(2)

x(i+1)=x(i)+h*(55*funkc(t(i),x(i))-59*funkc(t(i-1),x(i-1))+37*funkc(t(i-2),x(i-2))-9*funkc(t(i-3),x(i-3)))/24;

t(i+1)=t(i)+h;

i=i+1;

end

5. Экстраполяционный метод Адамса 4 порядка.

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

function [t,x] = Adamsa3_interpol(funkc,interval,inital,h)

t(1)=interval(1);

t(2)=t(1)+h;

t(3)=t(2)+h;

[bt,bx]= eilera_koshi_difur(funkc,[interval(1) interval(1)+2*h],inital,h/10);

x(1)=inital;

x(2)=bx(11);

x(3)=bx(21);

funkc=inline(funkc);

i=3;

while t(i)<interval(2)

t(i+1)=t(i)+h;

x(i+1)=(x(i)+19*h*funkc(t(i),x(i))/24-5*h*funkc(t(i-1),x(i-1))/24+h*funkc(t(i-2),x(i-2))/24+3*h/(8*t(i+1)*(t(i+1)-1)))/(1+3*h/(4*(t(i+1)-1)));

i=i+1;

end

Размещено на Allbest.ru