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

5. Построить отдельно графики , , , а также трехмерный график движения точки в декартовой системе координат средствами Matlab.

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

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

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

В процессе решения задачи оптимизации обычно необходимо найти оптимальные значения некоторых параметров, определяющих данную задачу. При решении инженерных задач их принято называть проектными параметрами. В качестве проектных параметров могут быть, в частности, значения линейных размеров объекта, массы, температуры и т.п. число 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 Интерполяция многочленом Лагранжа

Пусть задана функция y=f(x). Часто нахождение значений этой функции может оказаться трудоемкой задачей. Например, - параметр в некоторой сложной задаче, после решения которой определяется значение f(x), или f(x) измеряется в дорогостоящем эксперименте. В этих случаях можно получить небольшую таблицу значений функции, но прямое нахождение ее значений при большом количестве значений аргумента нереально. В такой ситуации f(x) заменяется приближенной функцией , которая в определенном смысле близка к функции f(x). Близость обеспечивается введением в функцию свободных параметров и их соответствующим выбором.

Итак, известны значения функции f(x) в точках , . Потребуем, чтобы для некоторой функции , где - свободные параметры, выполнялись равенства:

. (7)

Если (7) рассматривать как систему для определения , то этот способ называется интерполяцией (Лагранжевой).

Если зависит от нелинейно, то интерполяция нелинейная, иначе интерполяция линейная. В случае линейной интерполяции можно записать

, (8)

где - система линейно-независимых функций.

Подставим (8) в (7). Относительно получаем линейную систему уравнений:

, . (9)

Для однозначной разрешимости системы должно быть .

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

(10)

Система функций, удовлетворяющая условию (10), называется чебышевской.

Наиболее простой и (для многих случаев) удобной является система функций , . Функция при этом представляет собой многочлен степени (интерполяционный многочлен) с коэффициентами .

Система уравнений (9) в этом случае имеет вид:

, (11)

Определителем этой системы является отличный от нуля определитель Вандермонда:

.

Отсюда следует, что интерполяционный многочлен существует и единственен.

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

Если y С[a,b]- многочлен степени n, то - искомый интерполяционный многочлен степени ,т.к. .

Так как при , то y C[a,b] делится на для любых , т.е.

.

Так как , то .

Таким образом,

. (12)

Такая форма записи интерполяционного многочлена называется многочленом Лагранжа и обозначается, как правило, .

Существуют и другие формы записи того же самого интерполяционного многочлена.

Если обозначить , то (6) можно записать в виде .

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 - Сравнительная оценка решения задачи различными решателями

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

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

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

2.1 Алгоритм метода барьерных функций

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

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

или .

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

или .

Шаг 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.

Приложение

интегрирование программа matlab задача

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

function var7

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.1; % начальное значение x

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

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

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

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

tspan = [0 T];

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

%стандартная функция 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;

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

[t2,Y2,te2,ye2] = RK2(@dvar,tspan,h0,Y0,0.01);

figure(1);

plot(t2,Y2(:,1),'--r'), grid on;

xlabel('Время');

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

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

hold off;

figure(2);

plot(t2,Y2(:,2),'--r'), grid on;

xlabel('Время');

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

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

hold off;

figure(3);

plot(t2,Y2(:,3),'--r'), grid on;

xlabel('Время');

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

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

hold off;

figure(4);

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

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

hold off;

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

fprintf('t=%14.7f Y(1)=%14.7f Y(2)=%14.7f Y(3)=%14.7f\n',T/2,interp1(t,Y(:,1),T/2,'linear'),...

interp1(t,Y(:,2),T/2,'linear'),interp1(t,Y(:,3),T/2,'linear'));

fprintf('t2=%14.7f Y2(1)=%14.7f Y2(2)=%14.7f Y2(3)=%14.7f\n',T/2,interp1(t2,Y2(:,1),T/2,'linear'),...

interp1(t2,Y2(:,2),T/2,'linear'),interp1(t2,Y2(:,3),T/2,'linear'));

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

videoobj = VideoWriter('var7','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 = lagranzh(tt,B,t);

end

F = zeros(3,1);

F(1) = Y(1)^2-sin(Y(2))-Y(3);

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

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

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

% Метод Рунге-Кутты с априорным выбором шага интегрирования

% ----

% y0 - начальное значение y

% x - интервал (массив независимой переменной), на котором ищется решение ОДУ

% h0 - шаг интегрирования

% f - функция расчета правых частей ОДУ

% e - заданная точность

jmax = 1000;

h = h0;

j = 1;

while j < jmax

%

N = length(y0);

H = max(x);

n = ceil((H-x(1))/h);

y = zeros(N,n+1);

x(1) = min(x);

for k=1:1:N

y(k,1) = y0(k);

end

for i=1:1:n

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

end

for i=2:1:n+1

k1 = h.*f(x(i-1),y(:,i-1));

k2 = h.*f(x(i-1)+h/2,y(:,i-1)+k1/2);

k3 = h.*f(x(i-1)+h/2,y(:,i-1)+k2/2);

k4 = h.*f(x(i),y(:,i-1)+k3);

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

y(:,i) = y(:,i-1) + dy;

end

%-

L(j) = n+1;

if j ~= 1

for i3=1:1:N

M22(i3) = max(abs(y(i3,:)));

end

M2 = max(M22);

for i1=1:1:N

for i2=1:1:L(j-1)

Y2(i1,i2) = y(i1,2*i2-1)/M2;

end

end

end

if j ~= 1

for i3=1:1:N

M(i3) = max(abs(Y2(i3,:) - Y1(i3,:)));

end

end

for i3=1:1:N

M11(i3) = max(abs(y(i3,:)));

end

M1 = max(M11);

for i1=1:1:N

for i2=1:1:L(j)

Y1(i1,i2) = y(i1,i2)/M1;

end

end

if (j > 1) && (max(M) <= e)

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

end

h = h/2;

j = j+1;

end

x = x';

y = y';

te = x(n+1);

ye = y(n+1,:);

function yy = lagranzh(x,y,xx)

% интерполяция многочленом Лагранжа

N = length(x);

yy=zeros(size(xx));

for k=1:N

t=ones(size(xx));

for j=[1:k-1, k+1:N]

t=t.*(xx-x(j))/(x(k)-x(j));

end

yy = yy + y(k)*t;

end

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

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),

%

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

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

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

while k < kmax

options = optimset('MaxFunEvals',5000,'MaxIter',5000);

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

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

if abs(PP(k)) <= e || abs(PP(k)) == Inf

X1 = x1(k,:),

F1 = myfun(X1),

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

else

k = k + 1;

r = r/C;

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

end

end

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

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

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

% разработанная функция

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

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

function f = myfun(x)

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

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

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

function [c,ceq] = mycon(x)

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

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

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

ceq = [];

%

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

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

function F = mysupp(x)

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

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

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

f = myfun(x);

P = myfine(x);

F = f - P;

function P = myfine(x)

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

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

global r;

[c,ceq] = mycon(x);

P = r*sum(1./c);

% P = r*sum(log(-c));

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