Algorytmy i metody numeryczne

Ćwiczenie Nr 3

ROZWIĄZYWANIE NIELINIOWYCH RÓWNAŃ ALGEBRAICZNYCH

Zadania

1. Przeanalizować i przetestować program rozwiązywania nieliniowego równania algebraicznego metodą bisekcji, metodą siecznych i metodą Newtona (Przykład 1).

2. Przeanalizować i przetestować program wyznaczania pierwiastków zespolonych wielomianu x^4+1=0 (Przykład 2).

3. Przeanalizować i przetestować program analizy układu elektrycznego z diodą tunelową (Przykład 3).

4. Przeanalizować i przetestować program rozwiązania układu równań za pomocą dwuwymiarowej metody Newtona (Przykład 4).

5. Przeanalizować i przetestować program rozwiązywania układu trzech równań nieliniowych poprzez sprowadzenie tego układu do jednego równania względem częstotliwości i rozwiązanie tego równania za pomocą procedury "fzero" (Zadanie 1).

6. Przeanalizować i przetestować program wyznaczania wartości prądów w układzie elektrycznym (Zadanie 2).

7. Równanie
ma pierwiastki
. Stosując metodę
połowienia, obliczyć dodatni pierwiastek tego równania zaczynając od przedziału
. Ile iteracji należy wykonać, aby obliczyć pierwiastek z dokładnością do czterech miejsc dziesiętnych? Jaki jest maksymalny błąd po tej liczbie iteracji?

8. Równanie
ma tylko jeden pierwiastek rzeczywisty. Znaleźć ten pierwiastek metodą połowienia.

9. Znaleźć pierwiastek równania
w przedziale
. Ile trzeba wykonać iteracji, aby metodą połowienia otrzymać przybliżoną wartość pierwiastka z błędem nie przekraczającym
?

10. Metodą połowienia znaleźć dodatni pierwiastek równania
z dokładnością­ do 10^-2.

11. Metodą połowienia znaleźć wartość zmiennej
, dla której przecinają się wykresy funkcji
i
. Wymagana dokładność obliczeń wynosi 10^-4 .

12. Korzystając z metody połowienia znaleźć najmniejszy dodatni pierwiastek równania:
a)
, b)
, c)
.

Odpowiedzi.

7. Trzynaście iteracji daje wynik
1,4141844. Maksymalny błąd w każdym kroku, określony jako różnica między dwiema ostatnimi iteracjami, jest równy
. Dla
13 mamy (1/2)¹³
0,0001222.

8.
= -1,3247421.

9. 14 iteracji.

10.

11. x = 0,6191.

12. a) x-= 1,13226773, b) x = 0,4450, c) x = 0,92102480.

Programy dla przykładów

%%%%%%%%%%%%%%%%

% PRZYKЈAD 1 %

%%%%%%%%%%%%%%%%

clear all

close all

warning off

clc

% ROZWIҐZYWANIE NIELINIOWEGO RУWNANIA ALGEBRAICZNEGO x*x-x-2=0

% METODҐ BISEKCJI, METODҐ SIECZNYCH I METODҐ NEWTONA

%===============================================================================

% Metoda bisekcji wg sieci dziaіaс z rys.5.1

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

x_dokl=2; % rozwi№zanie dokіadne

a=1.5; b=3; % lewa i prawa granica przedziaіu poszukiwania rozwi№zania

dokladnosc=4e-15; % wskaџnik poї№danej dokіadnoњci rozwi№zania

deltax=(b-a)/2; % oszacowanie bікdu bezwglкdnego pierwszego przybliїenia

iter=0; % pocz№tek pкtli iteracyjnej (metody bisekcji)

while deltax>dokladnosc

iter=iter+1;

x=(a+b)/2;

if (f51(a)*f51(x)>0) % odwoіanie do zdefiniowanej niїej funkcji f51

a=x;

else b=x;

end

deltax=b-a; % oszacowanie bікdu bezwglкdnego kolejnego przybliїenia

dx_bis(iter)=abs(x_dokl-x)/x_dokl; % moduі bікdu wzglкdnego kolejnego

end % koniec pкtli iteracyjnej (metody bisekcji) przybliїenia

% Parametry rysunkуw

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

figsize=get(0,'ScreenSize'); % identyfikacja wielkoњci ekranu

figx=figsize(3)/2-10; if figx>630, figx=630;end % szerokoњж rysunkуw

figy=figsize(4)/2-100; if figy>412, figy=412;end % wysokoњж rysunkуw

figpos=[5 5 figx figy]; % poіoїenie rysunku

f=figure(1);

set(f,'Pos',figpos(1,:));

hold off

semilogy(dx_bis,':x')

hold on

% Metoda siecznych wg wzoru (5.6)

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

ximinus1=3;

xi=1.5;

deltax=abs(f51(x));

dokladnosc=1e-15; % wskaџnik poї№danej dokіadnoњci rozwi№zania

iter=0; % pocz№tek pкtli iteracyjnej (metody siecznych)

while deltax>dokladnosc

iter=iter+1;

xiplus1=xi-(xi-ximinus1)/(f51(xi)-f51(ximinus1))*f51(xi);

deltax=abs(f51(xiplus1));

ximinus1=xi;

xi=xiplus1;

dx_siecz(iter)=abs(x_dokl-xi)/x_dokl; % moduі bікdu wzglкdnego kolejnego

end % koniec pкtli iteracyjnej (metody siecznych) przybliїenia

semilogy(dx_siecz,':*r')

% Metoda Newtona wg wzoru (5.4)

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

xi=3;

deltax=abs(f51(xi));

dokladnosc=1e-15; % wskaџnik poї№danej dokіadnoњci rozwi№zania

iter=0; % pocz№tek pкtli iteracyjnej (metody Newtona)

while deltax>dokladnosc

iter=iter+1;

xiplus1=xi-f51(xi)/fprim51(xi); % odwoіanie do zdefiniowanej

xi=xiplus1; % niїej funkcji fprim51

deltax=f51(xi);

dx_newton(iter)=abs(x_dokl-xi)/x_dokl; % moduі bікdu wzglкdnego kolejnego

end % koniec pкtli iteracyjnej (metody Newtona) przybliїenia

% Prezentacja wynikуw

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

semilogy(dx_newton,':ok')

l=legend('metoda bisekcji','metoda siecznych','metoda Newtona');

set(l,'FontSize',8);

title('Zaleїnoњж moduіu bікdu wzglкdnego od liczby iteracji');

xlabel('Numer iteracji')

ylabel('Moduі bікdu wzglкdnego rozwi№zania');

grid on

fprintf('\n LICZBA ITERACJI\n')

fprintf(['METODA 1 2 3 4 5',...

' 6\n'])

fprintf('BISEKCJI %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n',dx_bis(1:6))

fprintf('SIECZNYCH %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n',dx_siecz(1:6))

fprintf('NEWTONA %10.2e %10.2e %10.2e %10.2e %6i%11i\n',[dx_newton(1:4),0,0])

function y=f51(x)

% obliczanie wartoњci funkcji definiuj№cej rуwnanie w przykіadzie 5.1

y=x.*x-x-2; % mnoїenie z kropk№, aby funkcja dziaіaіa poprawnie dla wektora arg.

function y=fprim51(x)

% obliczanie wartoњci pochodnej funkcji definiuj№cej rуwnanie w przykіadzie 5.1

y=2*x-1;

%%%%%%%%%%%%%%%%

% PRZYKЈAD 2 %

%%%%%%%%%%%%%%%%

clear all

close all

clc

format long

% WYZNACZANIE PIERWIASTKУW ZESPOLONYCH WIELOMIANU x^4+1=0

%===============================================================================

% Parametry rysunkуw

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

figsize=get(0,'ScreenSize'); % identyfikacja wielkoњci ekranu

figx=figsize(3)/2-10; if figx>630, figx=630;end % szerokoњж rysunkуw

figy=figsize(4)/2-100; if figy>412, figy=412;end % wysokoњж rysunkуw

figpos=[5 5 figx figy]; % poіoїenie rysunku

% WARIANT 1: rozkіad wielomianu na trуjmiany kwadratowe metod№ Bairstowa

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

% Punkt startowy: u0=1, v0=2

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

uv0=[1;2]; % punkt startowy

dok=1; % wskaџnik poї№danej dokіadnoњci rozwi№zania

uvi=uv0;

iter=0;

uinf=sqrt(2);vinf=1; % rozwi№zanie dokіadne

while dok>1e-15 % pocz№tek pкtli iteracyjnej realizowanej wg wzoru (5.19)

uvip1=uvi-muvp52(uvi)*fuvp52(uvi); % odwoіanie do zdefiniowanych niїej funkcji

% muvp52 i fuvp52

dok=norm(uvip1-uvi); % norma rуїnicy miкdzy dwoma kolejnymi przybliїeniami

uvi=uvip1;

iter=iter+1;

delta(iter)=sqrt((uvi(1)-uinf)^2+(uvi(2)-vinf)^2);

end % koniec pкtli iteracyjnej

% Prezentacja wynikуw dla punktu startowego u0=1, v0=2

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

%delta(find(delta<eps))=eps; % zast№pienie zer liczb№ eps

fprintf('\nu=%17.15f',uvi(1)); % (ze wzglкdu na wykres logarytmiczny)

fprintf(' v=%17.15f\n',uvi(2));

f=figure(1);

set(f,'Pos',figpos(1,:));

hold off

semilogy(1:iter,delta,'o-k')

h=xlabel('\iti');

set(h,'FontName','Times','FontSize',12);

h=ylabel('\it\Delta_i ');

set(h,'Rotation',0,'FontName','Times','FontSize',12);

t=title(['Zaleїnoњж bікdu od liczby iteracji dla punktu startowego ',...

'\itu\rm_0=1, \itv\rm_0=2']);

set(t,'FontSize',8);

grid on

fprintf('\nAby kontynuowaж, naciњnij klawisz\n')

pause

% Punkt startowy: u0=-2, v0=2

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

clear delta

uv0=[-2;2];

dok=1;

uvi=uv0;

iter=0;

uinf=-sqrt(2);

vinf=1;

while dok>1e-15 % pocz№tek pкtli iteracyjnej realizowanej wg wzoru (5.19)

uvip1=uvi-muvp52(uvi)*fuvp52(uvi);% odwoіanie do zdefiniowanych niїej funkcji

% muvp52 i fuvp52

dok=norm(uvip1-uvi);% norma rуїnicy miкdzy dwoma kolejnymi przybliїeniami

uvi=uvip1;

iter=iter+1;

delta(iter)=sqrt((uvi(1)-uinf)^2+(uvi(2)-vinf)^2);

end % koniec pкtli iteracyjnej

% Prezentacja wynikуw dla punktu startowego u0=-2, v0=2

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

%delta(find(delta<eps))=eps;% zast№pienie zer liczb№ eps

fprintf('\nu=%17.15f',uvi(1)); % (ze wzglкdu na wykres logarytmiczny)

fprintf(' v=%17.15f\n',uvi(2));

figure(1)

hold off

semilogy(1:iter,delta,'o-k')

h=xlabel('\iti');

set(h,'FontName','Times','FontSize',12);

h=ylabel('\it\Delta_i ');

set(h,'Rotation',0,'FontName','Times','FontSize',12);

t=title(['Zaleїnoњж bікdu od liczby iteracji dla punktu startowego ',...

'\itu\rm_0=-2, \itv\rm_0=2']);

set(t,'FontSize',8);

grid on

%===============================================================================

% WARIANT 2: wyznaczanie wartoњci wіasnych macierzy stowarzyszonej

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

C=[0 0 0 -1 % macierz stowarzyszona wg wzoru (5.20)

1 0 0 0

0 1 0 0

0 0 1 0];

fprintf('\nWyznaczone wartoњci wіasne macierzy stowarzyszonej\n')

eig(C)

format short

function [cdprim]=muvp52(uvi)

% Wyznaczanie macierzy odwrotnej do macierzy pochodnych wg wzoru (5.19)

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

ui=uvi(1);

vi=uvi(2);

J=[-3*ui.^2+2*vi,2*ui;-2*ui*vi,2*vi-ui.^2];

cdprim=inv(J);

function [cd]=fuvp52(uvi)

% Wyznaczanie wartoњci funkcji definiuj№cych ukіad rуwnaс (5.18)

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

ui=uvi(1);

vi=uvi(2);

cd=[2*ui.*vi-ui.^3;vi.^2-ui.^2.*vi+1];

%%%%%%%%%%%%%%%%

% PRZYKЈAD 3 %

%%%%%%%%%%%%%%%%

clear all

close all

clc

% ANALIZA UKЈADU ELEKTRYCZNEGO Z DIODҐ TUNELOWҐ

%===============================================================================

% Parametry rysunkуw

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

figsize=get(0,'ScreenSize'); % identyfikacja wielkoњci ekranu

figx=figsize(3)/2-10; if figx>630, figx=630;end % szerokoњж rysunkуw

figy=figsize(4)/2-100; if figy>412, figy=412;end % wysokoњж rysunkуw

figpos=[5 5 figx figy]; % poіoїenie rysunku

% WARIANT 1: sprowadzenie ukіadu rуwnaс do jednego rуwnania wzglкdem napiкcia

% na diodzie i rozwi№zanie tego rуwnania za pomoc№ procedury "fzero"

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

ud=fzero('f53',0); % odwoіanie do zdefiniowanej niїej funkcji f53

k=1;

jg=1;

R1=10;

id=k*ud*(ud*ud/3-3*ud/2+2); % wyznaczenie pr№du id dla znalezionego rozwi№z. Ud

udr=ud;

idr=id;

fprintf('\nRozwi№zanie uzyskane za pomoc№ procedury "fzero":\n')

fprintf(' ud=%6.4f id=%6.4f\n\n',ud,id)

fprintf('\nAby kontynuowaж, naciњnij klawisz\n')

pause

% WARIANT 2: rozwi№zanie ukіadu rуwnaс za pomoc№ procedury "fsolve"

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

% punkt startowy (2, 0.5)

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

%r1=fsolve('u53',[2;0.5],optimset('fsolve')); % odwoіanie do zdefiniowanej niїej funkcji "u53"

r1=fsolve('u53',[2;0.5]); % odwoіanie do zdefiniowanej niїej funkcji "u53"

ud=r1(1);

id=r1(2);

fprintf('\nRozwi№zanie uzyskane za pomoc№ procedury "fsolve"\n')

fprintf(' dla punktu startowego [2;0.5]:\n')

fprintf(' ud=%6.4f id=%6.4f\n\n',ud,id)

% punkt startowy (0.4, 0.05)

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

r2=fsolve('u53',[0.4,0.05],optimset('fsolve'));

ud=r2(1);

id=r2(2);

fprintf('\n dla punktu startowego [0.4;0.05]:\n')

fprintf(' ud=%6.4f id=%6.4f\n\n',ud,id)

fprintf('\nAby kontynuowaж, naciњnij klawisz\n')

pause

% WARIANT 3: rozwi№zanie ukіadu rуwnaс metod№ graficzn№

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

f=figure(1);

set(f,'Pos',figpos(1,:));

hold off

ud=(0:0.1:3);

id=k*ud.*(ud.*ud/3-3*ud/2+2);

hold off

plot(ud,id)

id=jg-ud/R1;

hold on

plot(ud,id,'.r')

h=line([udr udr],[0 idr]);

set(h,'Color','g','LineStyle',':');

h=line([0 udr],[idr idr]);

set(h,'Color','g','LineStyle',':');

h=xlabel('\itu_d\rm [V]');

set(h,'FontName','Times','FontSize',12);

h=ylabel('\iti_d\rm [mA]');

set(h,'Rotation',0,'FontName','Times','FontSize',12);

text(0.5,1,'ig-ud/R1')

h=line([r2(1) r2(1)],[0 r2(2)]);

set(h,'Color','k','LineStyle',':');

text(1.2,0.5,'minimum lokalne')

annotation('arrow',[0.5 r2(1)/2.68],[0.42 r2(2)/3*2]);

title('Graficzne rozwi№zanie ukіadu rуwnaс')

function y=f53(Ud)

% Wyznaczanie wartoњci funkcji reprezentuj№cej rуwnanie

% z jedn№ niewiadom№, do ktуrego zostaі sprowadzony ukіad rуwnaс

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

k=1;ig=1;R1=10;

id=k*Ud*(Ud*Ud/3-3*Ud/2+2);

y=Ud/R1+id-ig;

function y=u53(x)

% Wyznaczanie wartoњci dwуch funkcji definiuj№cych ukіad rуwnaс

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

k=1;ig=1;R1=10; % parametry ukіadu

ud=x(1); id=x(2); % zmiana nazw zmiennych wejњciowych uіatwiaj№ca interpretacjк

y(1,1)=id-k*ud*(ud*ud/3-3*ud/2+2);

y(2,1)=ud/R1+id-ig;

%%%%%%%%%%%%%%%%

% PRZYKЈAD 4 %

%%%%%%%%%%%%%%%%

clear all

close all

clc

% ROZWIҐZANIE UKЈADU RУWNAС ZA POMOCҐ DWUWYMIAROWEJ METODY NEWTONA

%===============================================================================

% Parametry rysunkуw

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

figsize=get(0,'ScreenSize'); % identyfikacja wielkoњci ekranu

figx=figsize(3)/2-10; if figx>630, figx=630;end % szerokoњж rysunkуw

figy=figsize(4)/2-100; if figy>412, figy=412;end % wysokoњж rysunkуw

figpos=[5 5 figx figy]; % poіoїenie rysunku

% Obliczenia

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

k=1;jg=1;R1=10; % parametry rуwnania

dokladnosc=1e-15; % wskaџnik poї№danej dokіadnoњci rozwi№zania

xinf=[2.3876725061573505;0.761232749384265]; % rozwi№zanie dokіadne

x=[5;0]; % punkt startowy

delta_i=1;

iter=0; % pocz№tek pкtli iteracyjnej (dwuwymiarowa metoda Newtona)

while delta_i>dokladnosc

iter=iter+1;

f=u53(x); % odwoіanie do zdefiniowanej niїej funkcji "u53"

fprim=fprim53(x); % odwoіanie do zdefiniowanej niїej funkcji "fprim53"

d_x=fprim\(-f); % rozwi№zanie ukіadu rуwnaс 5.23

x=x+d_x;

delta_i(iter)=sqrt((x(1)-xinf(1))^2+(x(2)-xinf(2))^2);

end

% Prezentacja wynikуw

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

delta_i(find(delta_i<eps))=eps;

fprintf('\nRozwi№zanie dwuwymiarow№ metod№ Newtona:\n')

fprintf('Ud=%6.4f id=%6.4f\n\n',x(1),x(2))

f=figure(1);

set(f,'Pos',figpos(1,:));

hold off

semilogy(delta_i,'*:k')

title('Bі№d rozwi№zania dwuwymiarow№ metod№ Newtona');

h=xlabel('\iti');

set(h,'FontName','Times','FontSize',12);

h=ylabel('||\Delta\it_i\rm||');

set(h,'Rotation',0,'FontName','Times','FontSize',12);

axis([0,10,1e-16,1e1]);

function y=u53(x)

% Wyznaczanie wartoњci dwуch funkcji definiuj№cych ukіad rуwnaс

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

k=1;ig=1;R1=10; % parametry ukіadu

ud=x(1); id=x(2); % zmiana nazw zmiennych wejњciowych uіatwiaj№ca interpretacjк

y(1,1)=id-k*ud*(ud*ud/3-3*ud/2+2);

y(2,1)=ud/R1+id-ig;

function y=fprim53(x)

% Obliczanie wartoњci pochodnych cz№stkowych dwуch funkcji,

% definiuj№cych ukіad rуwnaс (5.3), wzglкdem zmiennych Ud i id

k=1;jg=1;R1=10;

ud=x(1); id=x(2); % zmiana nazw zmiennych wejњciowych uіatwiaj№ca interpretacjк

y(1,1)=-k*(ud*ud-3*ud+2); % pochodna pierwszego rуwnania wzglкdem Ud

y(1,2)=1; % pochodna pierwszego rуwnania wzglкdem id

y(2,1)=1/R1; % pochodna drugiego rуwnania wzglкdem Ud

y(2,2)=1; % pochodna drugiego rуwnania wzglкdem id

%%%%%%%%%%%%%%%

% ZADANIE 1 %

%%%%%%%%%%%%%%%

clear all

close all

clc

format long

% ROZWIҐZYWANIE UKЈADU TRZECH RУWNAС NIELINIOWYCH

% POPRZEZ SPROWADZENIE TEGO UKЈADU DO JEDNEGO RУWNANIA WZGLКDEM CZКSTOTLIWOЊCI

% I ROZWIҐZANIE TEGO RУWNANIA ZA POMOCҐ PROCEDURY "fzero"

%===============================================================================

for xs=8e9:1e9:5e10 % wyznaczanie i pokazywanie kolejnych rozwi№zaс;

% eliminacja rozwi№zaс faіszywych

[x,f]=fzero('zad51',xs); % odwoіanie do zdefiniowanej niїej funkcji 'zad51'

if abs(f)<1e-3

fprintf('\nCzкstotliwoњж drgaс wіasnych fn=%7.4f [GHz]',x*1e-9);

fprintf('\nWartoњж funkcji dla uzyskanego rozwi№zania %10.4e\n',f);

fprintf('\nAby kontynuowaж, naciњnij klawisz\n')

pause

end

end

function z=zad51(f)

% Funkcja definiuj№ca rуwnanie nieliniowe wzglкdem czкstotliwoњci

% rуwnowaїne rozwi№zywanemu ukіadowi trzech rуwnaс

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

c=3e+8;

h=3e-3;

a=0.025;

L=0.025;

h=0.003;

epsilonr=10;

u1=(2*pi*f).^2/c/c;

u2=3.83171^2/a/a;

k2=u1*epsilonr-u2;

k=k2.^(0.5);

k02=u1-u2;

k0=k02.^(0.5);

z=sin(k*h).*cos(k0*(L-h))./k+sin(k0*(L-h)).*cos(k*h)./k0;

%%%%%%%%%%%%%%%

% ZADANIE 2 %

%%%%%%%%%%%%%%%

clear all

close all

clc

global R3 % wartoњж R3 przekazywana do uz52 jako zmienna globalna

% WYZNACZANIE WARTOЊCI PRҐDУW W UKЈADZIE ELEKTRYCZNYM PRZEDSTAWIONYM NA RYS.5.8

%===============================================================================

% WARIANT A: wyznaczanie wartoњci pr№dуw dla R3=0.1

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

R3=0.1;

options=optimset('fsolve');

options=optimset(options,'MaxIter',100000,'MaxFunEvals',100000,'TolFun',1e-12);

i=fsolve('uz52',[4,42,0.5],options); % odwoіanie do zdefiniowanej niїej funkcji

i1=i(1); % "uz52"

id=i(2)-i(1);

fprintf('\nRozwi№zanie uzyskane dla R3=%3.1f kOhm:\n',R3)

fprintf('id=%5.2f mA i1=%6.4f mA\n\n',id,i1)

% WARIANT B: wyznaczanie wartoњci pr№dуw dla R3=0

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

R3=0;

options=optimset('fsolve');

options=optimset(options,'MaxIter',100000,'MaxFunEvals',100000,'TolFun',1e-12);

i=fsolve('uz52',[0.8,142,0.8],options);

i1=i(1);

id=i(2)-i(1);

fprintf('\nRozwi№zanie uzyskane dla R3=%3.1f kOhm:\n',R3)

fprintf('id=%5.2f mA i1=%6.4f mA\n\n',id,i1)

3

---