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

Рисунок 5 - Оценка качества решения разработанной программой

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

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

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

2.1 Постановка задачи минимизации

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

Требуется найти локальный минимум целевой функции на множестве X, т.е. такую точку , что

,

где .

2.2 Алгоритм метода множителей

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

Приложени

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

function var4

clear all;

global T;

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

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

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

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

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

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

tspan = [0 T];

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

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

[t,Y] = 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] = 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;

% равномерная сетка

H = 0.5;

t1 = min(t):H:max(t);

% решения, найденные с помощью стандартного решателя Matlab

Y1(:,1) = lagrange(t,Y(:,1),t1);

Y1(:,2) = lagrange(t,Y(:,2),t1);

Y1(:,3) = lagrange(t,Y(:,3),t1);

% решения, найденные с помощью разработанной программы

Y2(:,1) = lagrange(t3,Y3(:,1),t1);

Y2(:,2) = lagrange(t3,Y3(:,2),t1);

Y2(:,3) = lagrange(t3,Y3(:,3),t1);

% оценка качества решения

Qx = simpson((Y1(:,1)-Y2(:,1)).^2,H),

Qy = simpson((Y1(:,2)-Y2(:,2)).^2,H),

Qz = simpson((Y1(:,3)-Y2(:,3)).^2,H),

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

videoobj = VideoWriter('var4','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)

F = zeros(3,1);

F(1) = Y(1)^2*tan(Y(2))-Y(3);

F(2) = (Y(1)+Y(2))^2-4*Y(3);

F(3) = Y(1)+Y(2)-Y(3);

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(1:end-1)';

y = y(1:end-1,:);

te = x_cur;

ye = y_cur';

function [yy] = lagrange(x,y,xx)

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

% x - массив координат узлов

% y - массив значений интерполируемой функции

% xx - массив значений аргумента, для которых надо вычислить значения полинома

% yy - массив значений полинома в точках xx

% ymax = max(y);

% y = y/ymax;

% узнаем число узлов интерполяции (N=n+1)

N=length(x);

% создаем нулевой массив значений интерполяционного полинома

yy=zeros(size(xx));

% в цикле считаем сумму по узлам

for k=1:N

% вычисляем произведения, т.е. функции Psi_k

t=ones(size(xx));

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

A = xx-x(j);

B = x(k)-x(j);

t = t.*A/B;

end

% накапливаем сумму

yy = yy + y(k)*t;

end

function [S] = simpson(y,h)

% вычисление определенного интеграла методом Симпсона

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

% y - значения функции в узловых точках

% S - рассчитанное значение определенного интеграла

n = length(y);

m = n/2;

sum1=0;

for i=2:1:m

sum1=sum1+y(2*i);

end

sum2=0;

for i=2:1:m-1

sum2=sum2+y(2*i+1);

end

S = h/3*(y(1)+y(n)+4*sum1+2*sum2);

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

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

function z2

% Метод множителей

clear all;

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

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

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

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

lambda(:,1) = xlsread('z2.xls',1,'D2:D3');

mu(1) = xlsread('z2.xls',1,'E2');

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

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

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

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

[X,Fval] = fmincon(@myfun,x(1,1:3),[],[],[],[],[],[],@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(@(x) mysupp(x,lambda(:,k),mu(k),r(k)),x(k,:));

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

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

if abs(PP(k)) <= e

X1 = x1(k,:),

F1 = myfun(X1),

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

elseif abs(PP(k)) > e

k = k + 1;

r(k) = C*r(k-1);

[c,ceq] = mycon(x1(k-1,:));

lambda(:,k) = lambda(:,k-1)+r(k-1).*ceq';

mu(k) = max(0,(mu(k-1)+r(k-1)*0));

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 = x(1)^2+2*x(2)^2+3*x(3)^2-2*x(1)*x(2)-2*x(1)*x(3);

function [c,ceq] = mycon(x)

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

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

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

ceq(1) = x(1)+x(2)+x(3)-1;

ceq(2) = x(1)-x(2)-2;

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

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

c = [];

function L = mysupp(x,lambda,mu,r)

% модифицированная функция Лагранжа

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

% lambda - параметр

% mu - параметр

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

f = myfun(x);

[c,ceq] = mycon(x);

P = myfine(x,mu,r);

S = lambda(1)*ceq(1)+lambda(2)*ceq(2);

L = f + S + P;

function P = myfine(x,mu,r)

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

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

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

% mu - параметр;

[c,ceq] = mycon(x);

P = r/2*sum(ceq.^2);

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