Iterační algoritmy

Algoritmus má cyklus, ve kterém se při každém jeho průchodu zpřesňuje výsledek.
Konec algoritmu je definován pomocí rozdílu mezi výsledky dvou následných iterací, který se díky zpřesňování výsledku v každé iteraci ve své absolutní hodnotě zmenšuje.

Iterační algoritmus pro výpočet určité matematické konstanty se setavuje pomocí matematické definice nekonečné řady, kterou je konstanta vyjádřena.
Algoritmus má cyklus, který kvůli neznámému počtu iterací nebude typu for ale while.
Cyklus má proměnnou, která má v následujících příkladech identifikátor "i", který tvoří určitou aritmetickou posloupnost.
Od této posloupnosti se odvozují výrazy, kterými se počítají konstanty jednotlivých členů nekonečné řady.
V nekonečné řadě
https://en.wikipedia.org/wiki/Madhava_of_Sangamagrama#The_value_of_%CF%80_(pi)
https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80
potřebujeme posloupnost lichých čísel 1, -3, 5, -7, 9, ...
Pro střídání znaménka zavedeme proměnnou "znamenko".
Aritmetická posloupnost tvořená identifikátorem "i" bude mít počáteční hodnotu 1 a přírůstek 2.
Přírůstek řady, jejíž členy se v absolutní hodnotě zmenšují, je její jednotlivý člen.
Matematická konstanta, kterou chceme iteračním algorimem vypočítat, je součtem jednotlivých členů řady, pro který zde máme proměnnou "suma".

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 0.01
int main (void)
{
  int i, znamenko = -1, jmenovatel;
  double suma = 0, prirustek;
  znamenko = -1;
  suma = 0;
  i = 1;
  do {
    znamenko = -znamenko;
    jmenovatel = znamenko * i;
    prirustek = 4.0 / jmenovatel;
    suma += prirustek;
    printf("%d\t", jmenovatel);
    printf("%.*f\t", POCET_DES_MIST, prirustek);
    printf("%.*f\n", POCET_DES_MIST, suma);
    i += 2;
  } while (fabs(prirustek) > PRESNOST);
  return 0;
}

Řada pro arcsin na stránce
https://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Infinite_series
konverguje mnohem rychleji, a proto můžeme v rozumném čase dosáhnout přesnějšího výsledku tím, že konstanta PRESNOST bude blíž nule.
Funkce arcsin(1/2) má výsledek 30 stupňů, tedy pi/6, viz rovnostranný trojúhelník se stranou rovnou 2, který se rozpůlí na dva pravoúhlé trojúhelníky s odvěsnami o délce 1, odmocnina ze 3 a přeponou o délce 2, viz
https://www.ck12.org/book/ck-12-algebra-ii-with-trigonometry-concepts/section/13.7/
Argument funkce arcsin má hodnotu 0.5 a v následujícím příkladu má identifikátor "z".
Členy matematické definice nekonečné řady si rozepíšeme a hledáme výraz, kterým z předhozího členu spočítáme další člen.
Člen č. 1 je     z
Člen č. 2 je     zzz    /(2*3)
Člen č. 3 je   3*zzzzz  /(2*4*5)
Člen č. 4 je 3*5*zzzzzzz/(2*4*6*7)
Podíly sousedních členů jsou
2. člen / 1. člen = (zzz/(2*3))             / z                   = 1*1*z*z/(2*3)
3. člen / 2. člen = (3*zzzzz/(2*4*5))       / (zzz/(2*3))         = 3*3*z*z/(4*5)
4. člen / 3. člen = (3*5*zzzzzzz/(2*4*6*7)) / (3*zzzzz  /(2*4*5)) = 5*5*z*z/(6*7)
Potřebujeme tedy aritmetickou posloupnost celých čísel s přírůstkem 2.
Každý člen nekonečné řady bude násobkem předchozího členu.

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-15
int main (void)
{
  int i;
  double prirustek, suma, z;
  prirustek = suma = z = 0.5;
  i = 1;
  do {
    prirustek *= i * i * z * z / ((i + 1) * (i + 2));
    suma += prirustek;
    printf("%.*f\t", POCET_DES_MIST, prirustek);
    printf("%.*f\t", POCET_DES_MIST, 6 * suma);
    printf("%g\n", 6 * suma - M_PI);
    i += 2;
  } while (prirustek > PRESNOST);
  return 0;
}

Sousední členy Fibonacciho řady konvergují ke konstantě zvané zlatý řez.
https://en.wikipedia.org/wiki/Golden_ratio
https://en.wikipedia.org/wiki/Fibonacci_number#Limit_of_consecutive_quotients

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-15
int main (void)
{
  int f1 = 0, f2 = 1, f3;
  double phi1, phi2, phi = (1.0 + sqrt(5.0)) / 2.0;
  do {
    f3 = f2 + f1;
    phi1 = (double) f3 / f2;
    phi2 = (double) f2 / f3;
    printf("%d\t%.*f\t%.*f\n", f3, POCET_DES_MIST, phi1, POCET_DES_MIST, phi2);
    f1 = f2;
    f2 = f3;
  } while (fabs(phi1 - phi2 - 1.0) > PRESNOST);
  printf("%.*f\t%g\n", POCET_DES_MIST, phi, phi - phi1);
  return 0;
}

Počáteční členy řady nemusí být 0 a 1 ale libovolné. Pokud jsou moc velké vzhledem ke konstantě PRESNOST, tak může dojít k přetečení datového typu int.

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-15
#define MAX_ITERATIONS 100
int main (void)
{
  int f1 = 300, f2 = 100, f3, i = 0;
  double phi1, phi2, phi = (1.0 + sqrt(5.0)) / 2.0;
  do {
    f3 = f2 + f1;
    phi1 = (double) f3 / f2;
    phi2 = (double) f2 / f3;
    printf("%d\t%.*f\t%.*f\n", f3, POCET_DES_MIST, phi1, POCET_DES_MIST, phi2);
    f1 = f2;
    f2 = f3;
    i++;
  } while (fabs(phi1 - phi2 - 1.0) > PRESNOST && i < MAX_ITERATIONS);
  printf("%.*f\t%g\n", POCET_DES_MIST, phi, phi - phi1);
  return 0;
}

Eulerovo číslo, základ přirozených logaritmů
https://en.wikipedia.org/wiki/E_(mathematical_constant)
Pro přesnost bližší nule například
#define PRESNOST 1e-8
přeteče datový typ int a poslední člen řady daný faktoriálem čísla 13 není správný, viz
https://en.wikipedia.org/wiki/Factorial

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-7
int main (void)
{
  int i, faktorial;
  double moje_e, prirustek;
  i = faktorial = 1;
  do {
    prirustek = 1.0 / faktorial;
    moje_e += prirustek;
    faktorial *= i;
    printf("%d\t%d\t%.*f\t%g\n", i, faktorial, POCET_DES_MIST, moje_e, moje_e - M_E);
    i++;
  } while (prirustek > PRESNOST);
  return 0;
}

Následující program počítá číslo e nejdříve od nejvyšších členů řady a potom ten samý počet členů řady sčítá znovu od nejnižšího členu řady.
Výsledná přesnost čísla e je nepatrně vyšší (4e-16 oproti 6e-16), když se sčítá řada od nejnižších členů, viz kapitola "1.3.1 Poučení" mé přednášky z předmětu Počítače I o sčítání od nejmenších hodnot.
V tomto programu je také použit datový typ long long int, jehož tisk v operačním systému Windows ukazují snímky "Datový typ long long" přednášky z jazyka C na stránce http://efis.tul.cz/~dana.nejedlova/ v odstavci "Programování I a II".

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-15
int main (void)
{
  int i;
  long long int faktorial;
  double moje_e = 0.0, prirustek;
  i = faktorial = 1;
  do {
    prirustek = 1.0 / faktorial;
    moje_e += prirustek;
    faktorial *= i;
    printf("%d\t%I64d\t%.*f\t%g\n", i, faktorial, POCET_DES_MIST, moje_e, moje_e - M_E);
    i++;
  } while (prirustek > PRESNOST);
  printf("\n");
  moje_e = 2.0;
  do {
    i--;
    printf("%d\t%I64d\t%.*f\t%g\n", i, faktorial, POCET_DES_MIST, moje_e, moje_e - M_E);
    faktorial /= i;
    moje_e += prirustek;
    prirustek = 1.0 / faktorial;
  } while (i > 1);
  return 0;
}

Výpočet druhé odmocniny
https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
Pro přesnost bližší nule například
#define PRESNOST 1e-14
a argument 1000 nekonverguje algoritmus k přesnosti, která by jej ukončila.

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-13
int main()
{
  double argument, odmocnina;
  printf("Zadej cislo, ktere chces odmocnit: ");
  scanf("%lf", &argument);
  odmocnina = argument / 2.0;
  do {
    odmocnina = (odmocnina + argument / odmocnina) / 2.0;
    printf("%.*f\t%.*f\t%g\n", POCET_DES_MIST, odmocnina, POCET_DES_MIST, argument / odmocnina, odmocnina * odmocnina - argument);
  } while (fabs(odmocnina * odmocnina - argument) > PRESNOST);
  return 0;
}

Metoda tečen pro výpočet jednoho z bodů, kde funkce protíná osu x.
https://en.wikipedia.org/wiki/Newton%27s_method

Vzorové zadání
Jaký je výsledek rovnice 3=cos(x)/ln(x-5)
Zdrojový kód programu pro toto zadání je:

#include <stdio.h>
#include <math.h>
#define POCET_DES_MIST 20
#define PRESNOST 1e-15
int main (void)
{
  double x0, x1, zmena;
  x0 = 6.1;
  do {
    x1 = x0 - (cos(x0) / log(x0 - 5.0) - 3.0) * log(x0 - 5) * log(x0 - 5) / (-sin(x0) * log(x0 - 5) - cos(x0) / (x0 - 5.0));
    printf("%.*f\n", POCET_DES_MIST, x1);
    zmena = x0 - x1;
    x0 = x1;
  } while (fabs(zmena) > PRESNOST);
  printf("\n%.*f\n", POCET_DES_MIST, cos(x1) / log(x1 - 5.0));
  return 0;
}


Konkrétně u tohoto příkladu platí, že počáteční hodnota 7, tedy, kdyby v programu bylo
x0 = 7.0;
způsobí, že se hodnota x0 změní na 4.145 a potom je argument logaritmu záporný a algoritmus selže.
V aplikaci MS Excel výchozí hodnota 7 nevadí, takže je zde použit trochu odlišný algoritmus, který však najde stejný kořen.
Pro zápis rovnice v jazyce C doporučuji dokumentaci knihovny math.h na stránce
http://www.cplusplus.com/reference/cmath/

Každý z vás dostane svoji úlohu odlišnou od ostatních a pokud chce započíst účast na jednom cvičení z Programování I, tak mi pošle:
1. Dokument s odvozením iteračního algoritmu pro svoji úlohu dle vzoru "6.docx" nebo "6.pdf".
2. Excelovský soubor, kde bude jeho úloha vypočtená pomocí nástroje Hledání řešení, viz soubor "6.xlsx". Tato úloha byla součástí testu z Excelu v zimním semestru předmětu Počítače I.
3. Zdrojový kód programu, který vypočte neznámou x v úloze, kterou jsem mu poslala, a hodnotu druhé strany rovnice, která se musí blížit zadané hodnotě (ve vzorovém zadání je to hodnota 3).