Page 9 - 6736
P. 9

%     N - максимальне число ітерацій
               %     eps - граничне відхилення f(x) у нулі
               %------------------------------------------------------------------------
               N=50;
               epselon=1e-7;
               eps=1e-7;
               %------------------------------------------------------------------------
               %Вихід:с - корінь рівяння f(x)=0
               %      Yc - значення функції f(x) при x=c
               %      err - похибка у знаходженні кореня
               %====================================
               %Функцію f(x) задати у підпрограмі fun_Solve
               %Вхід: х0 - початкова точка інтервалу
               %      xk - кінцева точка інтервалу
               %      delta -приріст аргументу функції f(x)
               %Вихід: Графік функції f(x)
               %====================================
               N=100;
               x0=0.51;
               xk=3;
               delta=(x0-xk)/N;
               %------------------------------------------------------------------------
               X=x0:delta:xk;
               q=2*c*(sqrt(1-(b./X).^2)+sqrt(1-m^2*(b./X).^2))./(pi*X);
               Y=X.*sqrt(q/(k*a1).*(a0*sqrt(k*q/a1-a2));
               plot(X,Y,'k-')
               %====================================
               %Із графіка функц. знаходимо інтервал
               %локального нул. Їх знач. вводимо
               %у командному вікні
               %====================================
               disp('Введіть а');
               a=input('a=');
               disp('Введіть b');
               b=input('b=');
               %------------------------------------------------
               Ya=fun_Solve(a,B,S);
               Yb=fun_Solve(b,B,S);
               if Ya*Yb>0
                   error('Неправильний інтервал локального нул.');
                   break
               end
               %====================================
               %Реаліз. процедури методу хорд
               %====================================
               for k=1:N
                   Dx=Yb*(b-a)/(Yb-Ya);
                   c=b-Dx;
                   ac=c-a;
                   Yc=fun_Solve(c,B,S);
                   if Yc==0,break
                   elseif Yb*Yc>0
                       b=c;
                       Yb=Yc;
                   else
                       a=c;
                       Ya=Yc;
                   end
                   Dx=min(abs(Dx),ac);
                   if abs(Dx)<epselon,break,end
                   if abs(Yc)<eps,break,end
               end
               omega=sqrt(a2/a0);
               %xsol=(b+a)/2;
               Yc=fun_Solve(c,B,S);
                                                               9
   4   5   6   7   8   9   10   11   12   13   14