Численные методы решения уравнений, + визуализация Опции

[junk]

Метод дихотомии (деления отрезка пополам)

Программа находит корень уравнения F(x) = 0, где F(x) - непрерывная на отрезке [a,b] функция, удовлетворяющая условию F(a)*F(B)<0. Для нахождения корня отрезок [a,b] делится пополам и выбирается тот полуинтервал, на концах которого знаки F(x) разные. Затем процесс деления повторяется до тех пор, пока длина интервала не станет меньше Eps.

{$n+}
Uses Crt;

Function f(x: Double): Double;
  Begin
    f := 1 / (Exp(x * Ln(2))) - 10 + 0.5*Sqr(x)
  End;

Var
  x, Eps, a, b, c: Double;
  n: Integer;

begin
  ClrScr;
  Writeln('Введите значения a и b'); Read(a, b);
  WriteLn('Введите точность Eps'); Read(Eps);

  n := 0;
  Repeat
    c := (a + b) / 2;
    If (f(a) * f©) < 0 Then b := c
    Else a := c;
    Inc(n)
  Until (b - a) <= Eps;

  x := (a + b) / 2;
  WriteLn('Корень равен  x=', x:10:7);
  WriteLn('Количество делений = ',n);
  ReadKey
end.

Метод хорд

{$n+}
Uses Crt;
Var
  x, Eps, a, b, c: Double;
  n: Integer;

Function f(x: Double): Double;
  Begin
    f := 1 / (Exp(x * Ln(2))) - 10 + 0.5*Sqr(x)
  End;

Function Chord(a, b: Double): Double;
  Begin
    Chord := (f(b)*a - f(a)*b) / (f(b) - f(a))
  End;

begin
  ClrScr;
  Writeln('Введите значения a и b');
  Read(a, b);

  WriteLn('Введите точность Eps');
  Read(Eps);

  ClrScr;
  n := 0;
  Repeat
    C := Chord(a, b);
    If f(a)*f© > 0 Then a := c
    Else b := c;
    Inc(n)
  Until Abs(Chord(a, b) - C) < Eps;

  x := c;
  WriteLn('Корень равен x=', x:10:7);
  WriteLn('Количество итераций =  ',n);
  ReadKey
end.

Метод Ньютона (касательных)

Действительный корень x уравнения F(x) = 0 вычисляется методом Ньютона по итерационному уравнению:
Xk+1 = Xk - F(Xk)/F'(Xk)

{$n+}
Function newton(start, Eps: Extended): Extended;
  Var
    X, prev: Extended;

  { Заданная функция }
  Function F(Arg: Extended): Extended;
    Begin
      F := Sqr(Arg) - 2
    End;

  { Производная заданной функции }
  Function Deriv(Arg: Extended): Extended;
    Begin
      Deriv := 2 * Arg
    End;

  Begin
    X := start;
    Repeat
      prev := X;
      X := prev - F(prev) / Deriv(prev);
    Until Abs(X - prev) <= Eps;
    newton := X
  End;

Var a, Eps: Extended;

begin
  WriteLn('Введите начальное приближение a');
  Read(a);
  WriteLn('Введите точность Eps');
  Read(Eps);

  WriteLn('Корень равен x= ', newton(a, Eps):10:7);
end.

Метод Ньютона с аппроксимацией производной

Если вычисление производной в методе Ньютона затруднено, можно заменить ее вычисление оценкой:
F'(X)= (F(X+h)-F(X))/h

{$n+}
Function newton(start, Eps: Extended): Extended;
  Var
    X, prev: Extended;

  Function noDeriv(Arg: Extended): Extended;

    Function F(Arg: Extended): Extended;
      Begin
        F := 1 / (Exp(arg * Ln(2))) - 10 + 0.5*Sqr(arg)
      End;

    Begin
      noDeriv := -2 * Eps * F(Arg) / (F(Arg + eps) - F(Arg - eps))
    End;

  Begin
    X := start;
    Repeat
      prev := X;
      X := prev + noDeriv(prev)
    Until Abs(X - prev) <= Eps;
    newton := X
  End;

Var a, Eps: Extended;

begin
  WriteLn('Введите начальное приближение a');
  Read(a);
  WriteLn('Введите точность Eps');
  Read(Eps);

  WriteLn('Корень равен x= ', newton(a, Eps):10:7);
end.

[SHnur]

Программа демонстрирует графический способ нахождения корней методом хорд!

Исходный код

program Metod_hord_graph;
uses
crt, graph;

var
gd, gm : integer;
maxY, minY, x1, x2, x3, y1, y2: real;

const
k3 = 3;{ Koeficenti }
k2 = -7;
k1 = -6;
k0 = 2; {}
grmaxX = 5; { Granici }
grminX = -5;

{---}
function f(x: real): real;
begin
f := k3*(x*x*x)+k2*(x*x)+k1*(x)+k0;
end;
{---}

function RealToScr(p: real; xory: boolean): integer;
var
i: integer;
ch: boolean;
rr: real;
begin
ch := false;
if xory then begin
rr := grminx;
for i := 0 to 645 do begin
if abs(rr-p) < 0.1 then begin
ch := true;
RealToScr := i;
end;
rr := rr + ((abs(grmaxx) + abs(grminx))/640 );
end;
end
else begin
rr := maxy;
for i := 0 to 479 do begin
if abs(rr-p) < 1 then begin
ch := true;
RealToScr := i;
end;
rr := rr - ((abs(maxy) + abs(miny))/480 );
end;
end;
if not(ch) then realtoscr := 10000;
end;
{---}

var
i, tx, ty: integer;
x: real;
s: string;

begin
{Graphic initialization}
gd := vga; gm := vgahi;
initgraph(gd,gm,'c:\programs\language\bp\bgi');
{}
maxY := f(grmaxX);
minY := f(grminX);
x1 := 1;
x2 := 5;
{ }

x := grminx;

for i := 0 to 639 do begin
putpixel(i,realtoscr(0,false),red);
end;

moveto(realtoscr(x,true),realtoscr(f(x),false));
repeat
x := x + 0.01;
tx :=realtoscr(x,true);
ty :=realtoscr(f(x),false);
if (tx <> 10000) and (ty <> 10000) then lineto(tx,ty);
delay(100);
until abs(x - grmaxx )<1E-5;

outtextxy(0,3,'Press any key !');
readkey;
setfillstyle(1,black);
bar(0, 0, 620, 10);

repeat
setcolor(green);
line(realtoscr(x1,true),realtoscr(f(x1),false),
realtoscr(x2,true),realtoscr(f(x2),false));
y1 := f(x1);
y2 := f(x2);
x3 := ((x1 * y2) - (x2 * y1)) / (y2 - y1);
if f(x3) < 0
then x1 := x3
else x2 := x3;
delay(10000);

until (abs(f(x3)) < 1);
str(x3 :2:5,s);
s := 'x = ' + s;
setcolor(white);
outtextxy(0,5,s);
readkey;
end.