Решение инженерных задач методами вычислительной математики

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

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ПЕНЗЕНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

Кафедра «Автоматика и телемеханика»

Курсовая работа

по дисциплине «Математические методы решения инженерных задач»

на тему «Решение инженерных задач методами вычислительной математики»

Выполнил студент: Джалолов С.Н

Группа: 15ПА2

Руководитель: к.т.н., доцент Козлов А.Ю.

2017

Содержание

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

Введение

1. Численное решение системы дифференциальных уравнений

2. Условная минимизация функций нескольких переменных

Заключение

Список использованных источников

Приложения

Введение

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

В курсовой работе рассмотрен метод Рунге-Кутты решения ОДУ.

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

В процессе решения задачи оптимизации обычно необходимо найти оптимальные значения некоторых параметров, определяющих данную задачу. При решении инженерных задач их принято называть проектными параметрами. В качестве проектных параметров могут быть, в частности, значения линейных размеров объекта, массы, температуры и т.п. число n проектных параметров x₁,x₂,…,x_n характеризует размерность (и степень сложности) задачи оптимизации.

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

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

1. Численное решение системы дифференциальных уравнений

1.1 Метод Рунге-Кутты

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

Конкретная прикладная задача может приводить к дифференциальному уравнению любого порядка или системе уравнений любого порядка. Однако обыкновенные дифференциальные уравнения n-го порядка можно с помощью замены свести к эквивалентной системе n уравнений первого порядка.

Различают три основных типа задач для решения обыкновенных дифференциальных уравнений: задачи Коши, краевые задачи и задачи на собственные значения. В этой курсовой работе рассмотрим метод решения задач Коши.

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

(1)

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

. (2)

Задача Коши состоит в том, чтобы найти функцию y=y(x), являющуюся решением уравнения (1) и удовлетворяющую условию (2).

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

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

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

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

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

(3)

Где

(4)

Заметим, что при использовании этого метода функцию необходимо вычислять четыре раза.

Отличительной чертой методов Рунге-Кутты является то, что при вычислении следующей точки x_m+1,y_m+1 используется информация только о точке x_m,y_m, но не предыдущих. В методах второго порядка и выше приходится вычислять значение функции в одной или нескольких промежуточных точках.

1.2 Выбор шага интегрирования

При численном решении задач Коши для ОДУ и систем ОДУ шаг численного решения можно выбирать апостериорно и априорно. В обоих случаях первоначальное значение шага h задается.

При апостериорном выборе шага последний изменяется в процессе счета на основе получаемой информации о поведении решения и на основе заданной точности .

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

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

2. Тем же методом с шагом решается задачи Коши с получением и .

3. Анализируется неравенство

. (5)

Если неравенство (5) удовлетворяется, то значение шага численного интегрирования на следующем шаге увеличивается вдвое по сравнению с первоначально выбранным шагом, т.е. становится равным 2h, и алгоритм повторяется, начиная с п.1.

4. Если неравенство (5) не выполняется, то счет ведется с шагом , начиная с отрезка и после получения значения анализируется неравенство

.

Если оно удовлетворяется, то дальнейший счет ведется с шагом и т.д.

При априорном выборе шага расчет ведется с первоначально выбранным шагом h с получением функции , , и с шагом с получением функции , Затем анализируется неравенство

. (6)

Если оно выполнено, то решение , , принимается за истинное, в противном случае расчет повторяется с шагом и сравниваются по норме (6) функции и и т.д.

1.3 Интерполяция многочленом Ньютона с разделенными разностями

Интерполяционная формула Ньютона позволяет выразить интерполяционный многочлен через значение в одном из узлов и через разделенные разности функции , построенные по узлам Она является разностным аналогом формулы Тейлора

(7)

Сначала приведем необходимые сведения о разделенных разностях. Пусть в узлах известны значения функции . Предположим, что среди точек нет совпадающих. Разделенными разностями первого порядка называются отношения

(8)

Будем рассматривать разделенные разности, составленные по соседним узлам, т. е. выражения

По этим разделенным разностям первого порядка можно построить разделенные разности второго порядка:

(9)

Аналогично определяются разделенные разности более высокого порядка. Например, если известны разности го порядка

то разделенная разность (+1)-го порядка определяется как

. (10)

Разделенная разность го порядка следующим образом выражается через значения функции в узлах:

(11)

Интерполяционным, многочленом Ньютона называется многочлен

(12)

1.4 Решение задачи

Для решения задачи интерполяции зададим в Excel исходные данные (рис. 1).

Рис. 1 Исходные данные для интерполяции

Для решения задачи интерполяции составили программу в Matlab. Текст программы приведен в Приложении А.1.

Для решения системы ОДУ методом Рунге-Кутты в Matlab применим функцию [t,Y,te,ye,ie]=ode45(@div,tspan,Y0,options) [2].

Текст программы приведен в Приложении А.1.

Результат решения системы ОДУ представлен на рисунках 2…4.

Рисунок 2 Изменение переменной x

Рисунок 3 Изменение переменной y



Рисунок 4 Изменение переменной z

Результат решения системы ОДУ стандартным и собственным решателем в декартовой системе координат представлен на рисунке 5.

Рисунок 5 Движение точки

В таблице 1 представлена сравнительная оценка решения задачи различными решателями.

Таблица 1

Сравнительная оценка решения задачи различными решателями

	Стандартный решатель	Разработанный решатель	, %
x(2.5)	11,4247208	11,4004557	0,2
y(2.5)	28,1827142	28,1462732	0,1
z(2.5)	8,1332118	8,1727713	0,5

Как видно из таблицы, погрешность решения задачи разработанным решателем не превышает 0,5%.

Для создания видеофайла решения задачи применим функцию videoobj = VideoWriter('v6','MPEG-4').

2. Условная минимизация функций нескольких переменных

2.1 Алгоритм метода штрафов

Шаг 1. Задать начальную точку ; начальное значение параметра штрафа ; число для увеличения параметра; малое число для остановки алгоритма. Положить .

Шаг 2. Составить вспомогательную функцию

.

Шаг 3. Найти точку безусловного минимума функции по с помощью какого-либо метода (нулевого, первого или второго порядка):

.

При этом задать все требуемые выбранным методом параметры. В качестве начальной точки взять . Вычислить .

Шаг 4. Проверить условие окончания:

а) если , процесс поиска закончить:

, ;

б) если , положить: , , и перейти к шагу 2.

2.2 Решение задачи

Составим и запустим программу в системе Matlab, реализующую метод штрафов для решения задачи условной минимизации функции нескольких переменных.

Текст программы приведен в Приложении А.2.

Результат выполнения программы представлен на рисунке 6.

Рисунок 6 Результат выполнения программы

Ввод исходных данных организован с помощью файла z2.xls (рис. 7).

Рисунок 7 Задание исходных данных в файле z2.xls

Вывод данных организован в файл z2.xls (рис. 8).

Рисунок 8 Вывод данных в файл z2.xls

Заключение

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

Список использованных источников

1. Дьяконов, В. Математические пакеты расширения Matlab. Специальный справочник [Текст] / В. Дьяконов, В. Круглов. Спб.: Питер, 2001. 480 с. 5000 экз. ISBN 5-318-0004-5.

2. Кетков, Ю.Л. MATLAB 7: программирование, численные методы [Текст] / Ю.Л. Кетков, А.Ю. Кетков, М.М. Шульц. Спб.: БХВ-Петербург, 2005. 752 с. 3000 экз. ISBN 5-94157-347-2.

3. Формалев, В.Ф. Численные методы [Текст] / В.Ф. Формалев, Д.Л. Ревизников. М.: ФИЗМАТЛИТ, 2004. 400 с. 1000 экз. ISBN 5-9291-0479-9.

4. Пантелеев, А.В. Методы оптимизации в примерах и задачах: Учебное пособие [Текст] / А.В. Пантелеев, Т.А. Летова. 2-е изд., исправл. М.: Высшая школа, 2005. 544 с. 3000 экз. ISBN 5-06-004137-9.

Приложение А.1

Текст программы решения задачи Коши для системы ОДУ

function var6

clear all;

global T;

global tt;

global B;

tt = xlsread('B.xls',1,'A2:A7');

B = xlsread('B.xls',1,'B2:B7');

T = 5; % время решения задачи

Y0(1) = 0.2; % начальное значение x

Y0(2) = 0; % начальное значение y

Y0(3) = 0.2; % начальное значение z

h0 = 0.1; % начальный шаг интегрирования

%--------------------------------

tspan = [0 T];

options = odeset('Events',@evar,'OutputFcn',@odeplot,'OutputSel',1,'Refine',4);

teout = [];

yeout = [];

ieout = [];

%стандартная функция Matlab------------------------------------------------

[t,Y,te,ye,ie] = ode45(@dvar,tspan,Y0,options);

odeplot([],[],'done'),grid;

figure(1);

plot(t,Y(:,1)), grid on; hold on;

xlabel('Время');

ylabel('переменная x');

figure(2);

plot(t,Y(:,2)), grid on; hold on;

xlabel('Время');

ylabel('переменная y');

figure(3);

plot(t,Y(:,3)), grid on; hold on;

xlabel('Время');

ylabel('переменная z');

figure(4);

plot3(Y(:,1),Y(:,2),Y(:,3)), grid on; hold on;

% собственная функция с апостериорным выбором шага-------------------------

[t3,Y3,te3,ye3] = RK1(@dvar,tspan,h0,Y0,0.001);

figure(1);

plot(t3,Y3(:,1),'--r'), grid on;

xlabel('Время');

ylabel('переменная x');

legend('Стандартный решатель','Разработанный решатель');

hold off;

figure(2);

plot(t3,Y3(:,2),'--r'), grid on;

xlabel('Время');

ylabel('переменная y');

legend('Стандартный решатель','Разработанный решатель');

hold off;

figure(3);

plot(t3,Y3(:,3),'--r'), grid on;

xlabel('Время');

ylabel('переменная z');

legend('Стандартный решатель','Разработанный решатель');

hold off;

figure(4);

plot3(Y3(:,1),Y3(:,2),Y3(:,3),'--r'), grid on;

legend('Стандартный решатель','Разработанный решатель');

hold off;

for k=1:1:length(t)

if t(k) >= T/2

fprintf('k=%d t=%14.7f Y(1)=%14.7f Y(2)=%14.7f Y(3)=%14.7f\n',k,t(k),Y(k,1),Y(k,2),Y(k,3));

break

end

end

for k=1:1:length(t3)

if t3(k) >= T/2

fprintf('k=%d t3=%14.7f Y3(1)=%14.7f Y3(2)=%14.7f Y3(3)=%14.7f\n',k,t3(k),Y3(k,1),Y3(k,2),Y3(k,3));

break

end

end

%--------------------------------------------------------------------------

teout = [teout; te];

yeout = [yeout; ye];

ieout = [ieout; ie];

%создание видеофайла-------------------------------------------------------

videoobj = VideoWriter('var6','MPEG-4');

open(videoobj);

fig = figure(5);

axis([min(Y(:,1)) max(Y(:,1)) min(Y(:,2)) max(Y(:,2)) min(Y(:,3)) max(Y(:,3))]);

axis manual;

hold on;

for k=1:1:length(t)

plot3(Y(k,1),Y(k,2),Y(k,3),'o','MarkerEdgeColor','g','MarkerFaceColor','g','MarkerSize',2,'LineWidth',2),grid on;

if k ~= 1

line([Y(k,1), Y(k-1,1)],[Y(k,2),Y(k-1,2)],[Y(k,3),Y(k-1,3)],'Color','g','LineStyle','-','LineWidth',4);% соединение точек линиями

end

F = getframe(fig);

writeVideo(videoobj,F);

end

close(fig);

close(videoobj);

%--------------------------------------------------------------------------

function [value,isterminal,direction] = evar(t,Y)

global T;

value = T - t;

isterminal = 1;

direction = -1;

function F = dvar(t,Y)

global tt;

global B;

if t <= tt(1)

Bt = B(1);

elseif t >= tt(end)

Bt = B(end);

else

Bt = nuton(tt,B,t);

end

F = zeros(3,1);

F(1) = sin(Y(1))*cos(Y(2))+Y(3);

F(2) = abs((Y(1)^2+Y(2)^2)*Y(2))^(1/8)+2*Y(3);

F(3) = Y(3)^2*cos(Y(2))+Bt;

function [ye] = RK(f, h, x_c, y_c)

k1 = h.*f(x_c,y_c);

k2 = h.*f(x_c + h/2, y_c + k1./2);

k3 = h.*f(x_c + h/2, y_c + k2./2);

k4 = h.*f(x_c + h, y_c + k3);

ye = y_c + (k1 + 2.*k2 + 2.*k3 + k4)./6;

function [t,y,te,ye] = RK1(f,x,h0,y0,e)

h = h0;

t(1) = x(1);

y(1,:) = y0;

k = 1;

flag = 0;

x_cur = x(1);

y_cur = y0';

while x_cur <= max(x)

if flag > 1

break;

end

[ye1] = RK(f, h, x_cur, y_cur);

[ye2] = RK(f, h/2, x_cur, y_cur);

[ye3] = RK(f, h/2, x_cur + h/2, ye2);

if max(abs(ye1 - ye3)) > e

h = h/2;

continue;

end

[ye11] = RK(f, h, x_cur + h, ye1);

[ye4] = RK(f, 2*h, x_cur, y_cur);

if max(abs(ye11 - ye4)) < e

h = h*2;

continue;

end

if (x_cur + h) > max(x)

h = abs(max(x) - x_cur);

end

x_cur = x_cur + h;

y_cur = ye1;

y(k+1,:) = y_cur;

t(k+1) = x_cur;

if x_cur >= max(x)

flag = flag + 1;

end

k = k + 1;

end

t = t';

te = x_cur;

ye = y_cur';

function yy = nuton(x,y,xx)

% интерполяция многочленом Ньютона

N = length(x);

DIFF = y;

for k = 1: N-1

for i = 1: N - k

DIFF(i) = (DIFF(i+1) - DIFF(i)) / (x(i+k) - x(i));

end

end

yy = DIFF(1) * ones(size(xx));

for k = 2: N

yy = DIFF(k) + (xx - x(k)).* yy;

end

Приложение А. 2

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

function z2

% Метод штрафов

clear all;

global r;

% исходные данные

x(1,:) = xlsread('z2.xls',1,'A2:A3');

r = xlsread('z2.xls',1,'B2');

C = xlsread('z2.xls',1,'C2');

e = xlsread('z2.xls',1,'D2');

%--------------------------------------------------------------------------

% минимизация стандартной функцией

options = optimoptions(@fmincon,'Algorithm','active-set','MaxIter',500,'Display','off');

[X,Fval] = fmincon(@myfun,x(1,:),[],[],[],[],[],[],@mycon,options),

xlswrite('z2.xls',Fval,2,'C2');

xlswrite('z2.xls',X',2,'C3');

%--------------------------------------------------------------------------

k = 1; % Счетчик шагов

kmax = 1000; % Предельное число шагов

while k < kmax

[x1(k,:),fval1(k)] = fminsearch(@mysupp,x(k,:)); % нахождение безусловного минимума вспомогательной функции

PP(k) = myfine(x1(k,:)); % расчет значения функции штрафа при оптимальных х на текущей итерации

if PP(k) <= e

X1 = x1(k,:),

F1 = myfun(X1),

break % Выход из цикла в случае выполнения условия

else

k = k + 1;

r = C*r;

x(k,:) = x1(k-1,:);

end

end

xlswrite('z2.xls',F1,2,'B2');

xlswrite('z2.xls',X1',2,'B3');

function f = myfun(x)

% целевая функция

% x - аргументы целевой функции

f = 4*x(1)^2+x(2)^2+2*x(2);

function [c,ceq] = mycon(x)

% расчет ограничений

% х - аргументы целевой функции

% ограничения типа "равенство"

ceq = [];

%--------------------------------------------------------------------------

% ограничения типа "неравенство"

c(1) = x(1)^2+x(2)^2+2*x(2);

c(2) = -x(1)+x(2)-2;

function F = mysupp(x)

% вспомогательная функция

% х - аргументы целевой функции

% r - значение параметра штрафа на текущей итерации

global r;

f = myfun(x);

P = myfine(x);

F = f + r/2*P;

function P = myfine(x)

% штрафная функция

% х - аргументы целевой функции

[c,ceq] = mycon(x);

P = sum(ceq.^2)+sum(max(0,c).^2);

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

...