Решение трансцендентного уравнения

Балтийский  государственный технический университет  «Военмех» им. Д.Ф. Устинова 
 
 
 
 
 

Лабораторная  работа №1

Решение трансцендентного уравнения

Цель  работы: Знакомство со средой программирования Matlab. Метод хорд и метод половинного деления. 
 
 
 
 
 
 
 
 
 

Выполнил: 

Преподаватель:

Овчинникова О.К. 
 
 
 
 
 
 
 
 
 
 

2010 г. 
 

% Главная программа  решения трансцендентного уравнения

% Контур сопла

x0 = 0; x1 = 1; x2 = 3;

y0 = 1.5; y1 = 1; y2 = 2;

N=50;

x = x0 : (x2 - x0)/(N - 1) : x2;

for k = 1:N

     if x(k) < x1

            y(k) = y0+(y1-y0)./(x1-x0).*(x(k)-x0);

      else

            y(k) = y1+(y2-y1)./(x2-x1).*(x(k)-x1);

      end

end

figure(1), clf

plot(x, y)

title(‘контур сопла’, FontName, ArialCyr, 14)

xlabel(‘x’)

ylabel(‘y’)

Fk = pi.*y1.^2;

F = pi.*y.^2;

Fk_F = Fk./F;

gam = 1.4; eps = [0.1  0.01  0.001];

for m = 1:3

      for k = 1:N

            if x(k) < x1

                  [t(k), M(k)] = solve1(Fk_F(k), gam, 0, eps(m));

            else

                  [t(k), M(k)] = solve1(Fk_F(k), gam, 1, eps(m));

            end

      end

end

t_iter(m) = mean(t);

figure(2), clf

plot(t_iter(m), eps)

title(‘Зависимость точности от количества итераций’, FontName, ArialCyr, 14)

xlabel(‘t’)

ylabel(‘eps’)

figure(3), clf

plot(M, q_fun(M, gam), ‘r’, x, M, ‘g’)

legend(‘q(M)’, ‘M(x)’)

figure (4), clf

plot(M, pi_fun(M, gam), ‘m’, M, eps_fun(M, gam), ‘y’, M, tau_fun(M, gam), ‘b’)

legend(‘pi’, ‘eps’, ‘tau’)

for k = 1:N

      if x(k) < x1

            M(k) = solve2(Fk_F(k), gam, 0 , eps(m));

      else

            M(k) = solve2(Fk_F(k), gam, 1, eps(m));

      end

end 

function q = q_fun(M, gam)

% Функция приведенного  расхода

% Аргументы М и  gam, возвращающие значения q_fun(M, gam)

q = M.*(2/(gam+1)*(1+(gam-1)/2.*M.^2)).^(-(gam)/2/(gam-1));

return 

function pi = pi_fun(M, gam)

% Первая функция  изоэнтропического течения

% Аргументы М и  gam

pi = (1+(gam-1)/2.*M.^2)).^((-gam)/(gam-1));

return 

function eps = eps_fun(M, gam)

% Вторая функция изоэнтропического течения

% Аргументы М и  gam

eps = (1+(gam-1)/2.*M.^2)).^((-1)/(gam-1));

return 

function tau = tau_fun(M, gam)

% Третья функция изоэнтропического течения

% Аргументы М и  gam

tau = (1+(gam-1)/2.*M.^2)).^(-1);

return 

function [M, t] = solve1(Fk_F, gam, key, eps, t)

% Метод дихотомии,  возвращает значение M

% Аргументы отношение Fk_F = Fk./F, gam, eps, t

% key = 1 – решение в сверхзвуке, key = 0 – решение в дозвуке

if key==1

      Ml = 1; Mr = 100;

else

      Ml = 0; Mr = 0.99;

end

ql = q_fun (Ml, gam)-Fk_F;

qr = q_fun (Mr, gam)-Fk_F;

t = 0;

while abs(Mr-Ml) > eps

      Mc = (Mr+Ml)/2;

      qc = q_fun(Mc, gam)-Fk_F;

      if qc*ql < 0

            Mr = Mc; qr = qc;

      else

            Ml = Mc; ql = qc;

      end

t = t+1

end

M = Mc;

return 

function M = solve2(Fk_F, gam, key, eps)

% Метод хорд, возвращает значение M

% Аргументы отношение Fk_F = Fk./F, gam, eps

% key = 1 – решение в сверхзвуке, key = 0 – решение в дозвуке

if key==1

      Ml = 1; Mr = 100;

else

      Ml = 0; Mr = 0.99;

end

ql = q_fun (Ml, gam)-Fk_F;

qr = q_fun (Mr, gam)-Fk_F;

while abs(Ml-Fk_F) > eps

      Mc = Ml-(ql/(ql-qr).*(Ml-Mr);

      qc = q_fun(Mc, gam)-Fk_F;

      if Mr > Ml

            Ml = Mc; ql = qc;

      else

            Mr = Mc; qr = qc;

      end

end

M = Mc;

return