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

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

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

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

 

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

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

 

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

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

 

Заключение

 

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

 

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

 

1. Дьяконов В., Круглов  В. Математические пакеты расширения  Matlab. Специальный справочник – Спб.: Питер, 2001, 480 с.

2. Кетков Ю.Л., Кетков А.Ю., Шульц М.М. MATLAB 7: программирование, численные методы - Спб.: БХВ-Петербург, 2005, 752 с.

3. Банди Б. Методы оптимизации. М. Радио и связь. 1988.

4. Формалев В.Ф., Ревизников Д.Л. Численные методы. – М.: ФИЗМАТЛИТ, 2004, 400 с.

5. Пантелеев А.В., Летова Т.А. Методы оптимизации в примерах и задачах. – М.: Высшая школа, 2005, 544 с.

 

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

 

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

function var24

clear all;

 

global T;

 

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

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

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

Y0(3) = 0.1; % начальное значение 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;

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

[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;

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

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('t3=%14.7f Y3(1)=%14.7f  Y3(2)=%14.7f  Y3(3)=%14.7f\n',T/2,interp1(t3,Y3(:,1),T/2,'linear'),...

        interp1(t3,Y3(:,2),T/2,'linear'),interp1(t3,Y3(:,3),T/2,'linear'));

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

videoobj = VideoWriter('var24','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)+sqrt(abs(Y(2)))-Y(3).^2;

F(2) =sqrt(abs(tan(Y(1)+Y(2))))+4*Y(3);

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

 

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';

 

 

Приложение А.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)+9/x(2)+x(1)+x(2);

 

function [c,ceq] = mycon(x)

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

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

 

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

ceq = [];

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

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

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

c(2) = -x(1)

c(2) = -x(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);