- #include<iostream>
- #include<iomanip>
- #include<conio.h>
- #include<cmath>
- #include<ctime>
- #define zloty = 0.618033988749
- double interpolacja_liniowa();
- double interpolacja_kwadratowa();
- double interpolacja_Lagrangea();
- double Cubic_Spline();
- double Wmaciez_CS(int size, double * v, double * u, double * w, int zmiana);//funkcja do licznia wyznacznika maciezy do Cubic_Spline() zmiana(0)-głowny wyznacznik,zmiana(1)-kolumny podmienione na v(i)
- void mciez_troj(int size, double ** maciez);
- //rozwiazywanie rownan
- double metoda_bisekcji(double(*F)(double));
- double metoda_zgadywana(double(*F)(double));//iteraqcyjna
- double zloty_podzial_min(double(*F)(double));
- double zloty_podzial(double(*F)(double));
- double metoda_falsi(double(*F)(double));
- double metoda_siecznych(double(*F)(double));
- double metoda_newtona_rapshona(double(*F)(double), double(*FP)(double));
- double metoda_newtona_rapshona_2_stopnia(double(*F)(double), double(*FP)(double), double(*FPP)(double));
- double metoda_halleya(double(*F)(double), double(*FP)(double), double(*FPP)(double));
- //funkcje matematyczne
- double funkcja(double x);
- double funkcja1(double x);
- double funkcja2(double x);
- double funkcja3(double x);
- double funkcja4(double x);
- double funkcja4P(double x);
- double funkcja4PP(double x);
- double funkcja5(double x);
- double funkcja6(double x);
- //calki
- double motoda_trapezow(double(*F)(double));
- double metoda_simpsona(double(*F)(double));
- double metoda_monte_carlo(double(*F)(double));//losowa
- int main(int argc, char* argv[]) {
- do {
- //std::cout << "interpolacja_liniow: " << interpolacja_liniowa();
- //std::cout << "ninterpolacja_kwadratowa: " << interpolacja_kwadratowa();
- //std::cout << "ninterpolacja_Lagrangea: " << interpolacja_Lagrangea();
- //std::cout << "nCubic_Spline: " << Cubic_Spline();
- //std::cout << "nnZloty Podzial min: " << zloty_podzial_min(funkcja2);
- //std::cout << "nnZloty Podzial: " << zloty_podzial(funkcja);
- //std::cout << "nnmetoda bisekcji: " << metoda_bisekcji(funkcja2);
- //std::cout << "nnmetoda zgadywana: " << metoda_zgadywana(funkcja1);
- //std::cout << "nnmetoda_falsi: " << metoda_falsi(funkcja4);
- //std::cout << "nnmetoda_siecznych: " << metoda_siecznych(funkcja4);
- //std::cout << "nnmetoda_newtona_rapshona: " << metoda_newtona_rapshona(funkcja4, funkcja4P);
- //std::cout << "nnmetoda_newtona_rapshona_2_stopnia: " << metoda_newtona_rapshona_2_stopnia(funkcja4, funkcja4P,funkcja4PP);
- //std::cout << "nnmetoda_halleya: " << metoda_halleya(funkcja4, funkcja4P, funkcja4PP);
- //std::cout << "nnmotoda_trapezow: " << motoda_trapezow(funkcja6);
- //std::cout << "nnmotoda_trapezow: " << metoda_simpsona(funkcja6);
- //std::cout << "nnmetoda_monte_carlo: " << metoda_monte_carlo(funkcja6);
- ..std::cout << "nn";
- std::cout << "end";
- if (getch() == 27)break;
- std::cout << "nn";
- } while (1);
- return 0;
- }
- double interpolacja_liniowa()
- {
- double x[3];
- double y[3];
- std::cout << "podaj punkt :n x: ";
- std::cin >> x[2];
- //std::cout << "y: ";
- //std::cin >> y[2];
- std::cout << "podaj punkt 0 :n x0: ";
- std::cin >> x[0];
- std::cout << "y0: ";
- std::cin >> y[0];
- std::cout << "podaj punkt 1 :n x1: ";
- std::cin >> x[1];
- std::cout << "y1: ";
- std::cin >> y[1];
- return y[2] = y[0] + ((y[1] - y[0]) / (x[1] - x[0]))*(x[2] - x[0]);
- }
- double interpolacja_kwadratowa()
- {
- double x[4];
- double y[4];
- std::cout << "podaj punkt :n x: ";
- std::cin >> x[3];
- //std::cout << "y: ";
- //std::cin >> y[3];
- std::cout << "podaj punkt 0 :n x0: ";
- std::cin >> x[0];
- std::cout << "y0: ";
- std::cin >> y[0];
- std::cout << "podaj punkt 1 :n x1: ";
- std::cin >> x[1];
- std::cout << "y1: ";
- std::cin >> y[1];
- std::cout << "podaj punkt 2 :n x2: ";
- std::cin >> x[2];
- std::cout << "y2: ";
- std::cin >> y[2];
- return y[3] = ((x[3] - x[1]) / (x[0] - x[1]))*((x[3] - x[2]) / (x[0] - x[2]))*y[0] + ((x[3] - x[0]) / (x[1] - x[0]))*((x[3] - x[2]) / (x[1] - x[2]))*y[1] + ((x[3] - x[0]) / (x[2] - x[0]))*((x[3] - x[1]) / (x[2] - x[1]))*y[2];
- }
- double interpolacja_Lagrangea()
- {
- int _size;
- std::cout << "ile :";
- std::cin >> _size;
- double *x = new double[_size + 1];
- double *y = new double[_size];
- double *s = new double[_size];
- std::cout << "podaj x do obliczen";
- std::cin >> x[_size];
- for (int i = 0; i < _size; i++) {
- std::cout << "podaj x" << i << ": ";
- std::cin >> x[i];
- std::cout << "podaj y" << i << ": ";
- std::cin >> y[i];
- }
- for (int i = 0; i < _size; i++)
- {
- s[i] = y[i];
- for (int j = 0; j < _size; j++)
- {
- if (x[i] != x[j])
- {
- s[i] *= (x[_size] - x[j]);
- }
- }
- for (int j = 0; j < _size; j++)
- {
- if (x[i] != x[j])
- {
- s[i] /= (x[i] - x[j]);
- }
- }
- }
- double suma = 0;
- for (int i = 0; i < _size; i++)
- {
- suma += s[i];
- }
- return suma;
- delete[] x,y,s;
- }
- double Cubic_Spline()
- {
- int _size = 0;
- std::cout << "ile :";
- std::cin >> _size;
- double mainx = 0;
- double mainy = 0;
- double W;
- std::cout << "dla jakiego x :";
- std::cin >> mainx;
- double *x = new double[_size];
- double *y = new double[_size];
- double *a = new double[_size];
- double *b = new double[_size];
- double *c = new double[_size];
- double *d = new double[_size];
- double *h = new double[_size];
- double *v = new double[_size];
- double *u = new double[_size];
- double *w = new double[_size];
- c[0] = 0;
- c[_size - 1] = 0;
- for (int i = 0;i < _size;i++)
- {
- std::cout << "podaj x" << i << " :";
- std::cin >> x[i];
- std::cout << "podaj y" << i << " :";
- std::cin >> y[i];
- }
- for (int i = 0; i < _size; i++)
- {
- a[i] = y[i];
- //std::cout << "a[" << i << "]" << a[i] << " ";
- }
- //std::cout << "nn";
- for (int i = 0; i < _size - 1; i++)
- {
- h[i] = x[i + 1] - x[i];
- //std::cout << "h[" << i << "]" << h[i] << " ";
- }
- //std::cout << "nn";
- for (int i = 1; i < _size - 1; i++)
- {
- u[i] = (h[i - 1] / (h[i - 1] + h[i]));
- //std::cout << "u[" << i << "]" << u[i] <<" ";
- }
- std::cout << "nn";
- for (int i = 1; i < _size - 1; i++)
- {
- w[i] = (h[i] / (h[i - 1] + h[i]));
- //std::cout << "w[" << i << "]" << w[i] << " ";
- }
- std::cout << "nn";
- for (int i = 1; i < _size - 1; ++i)
- {
- v[i] = (((y[i + 1] - y[i] / h[i]) - (y[i] - y[i - 1] / h[i - 1])) / (h[i - 1] + h[i]));
- //std::cout << "v[" << i << "]" << v[i] << " ";
- }
- W = Wmaciez_CS(_size,v,u,w,0);
- //std::cout << "nn wyznacznik = " << W << "nn";
- for (unsigned int i = 1; i < _size - 1; ++i)
- {
- c[i] = (Wmaciez_CS(_size, v, u, w, i) / W);
- //std::cout << "c[" << i << "]" << c[i] << " ";
- }
- //std::cout << "nn";
- //for (unsigned int i = 1; i < _size - 1; ++i)
- //{
- // std::cout << "w[" << i << "]" << Wmaciez_CS(_size, v, u, w, i) << " ";
- //}
- //std::cout << "nn";
- for (int i = 0; i < _size - 1; ++i)
- {
- d[i] = ((c[i + 1] - c[i]) / 3 * h[i]);
- //std::cout << "d[" << i << "]" << d[i] << " ";
- }
- //std::cout << "nn";
- for (unsigned int i = 0; i < _size; i++)
- {
- b[i] = (1 / h[i]) * (y[i + 1] - y[i] - c[i] * (h[i] * h[i]) - d[i] * (h[i] * h[i] * h[i]));
- if (i == _size - 1)
- {
- b[i] = b[i - 1] + 2 * (c[i - 1] * h[i - 1]) + 3 * (d[i - 1] * h[i - 1]);
- }
- //std::cout << "b[" << i << "]" << b[i] << " ";
- }
- std::cout << "nn";
- //x mniejsze
- if (mainx < x[0])
- {
- mainy = a[0] + b[0] * (mainx - x[0]);
- }
- //x wieksze
- if (mainx > x[_size-1])
- {
- mainy = a[_size - 1] + b[_size-2] * (mainx - x[_size - 1]);
- }
- //x w przedziale
- for (int i = 0; i < _size - 1; i++)
- {
- if (x[i] <= mainx && mainx <= x[i + 1])
- {
- mainy = a[i] + b[i] * (mainx - x[i]) + c[i] * ((mainx - x[i])*(mainx - x[i])) + d[i] * ((mainx - x[i])*(mainx - x[i])*(mainx - x[i]));
- }
- }
- delete[] x,y,a,b,c,d,w,v,u,h;
- return mainy;
- }
- double Wmaciez_CS(int size, double * v, double * u, double * w, int zmiana)
- {
- double **maciez = new double *[size];
- //wyznacznik
- double W = 1.0;
- size -= 2;
- for (unsigned int i = 0; i < size; i++)
- {
- maciez[i] = new double[size];
- }
- for (unsigned int i = 0; i < size; i++)
- {
- for (unsigned int j = 0; j < size; j++) {
- if (i == j) maciez[i][j] = 2;
- else if (i + 1 == j) maciez[i][j] = u[j];
- else if (i == j + 1) maciez[i][j] = w[j + 1];
- else maciez[i][j] = 0;
- }
- if (zmiana != 0)
- {
- for (unsigned int i = 0; i < size; i++)
- {
- maciez[i][zmiana - 1] = 3 * v[i + 1];
- }
- }
- }
- for (int w = 0; w < size; w++) {
- mciez_troj(size, maciez);//licze maciez trójkątną
- }
- for (int i = 0; i < size; i++)
- {
- W *= maciez[i][i];
- }
- for (unsigned int i = 0; i < size; i++)
- {
- delete[] maciez[i];
- }
- delete[] maciez;
- return W;
- }
- void mciez_troj(int size, double ** maciez)
- {
- for (int i = 0; i < size - 1; i++)
- {
- for (int j = i + 1; j < size; j++)
- {
- for (int q = size - 1; q >= i; q--)
- {
- maciez[j][q] = maciez[j][q] - maciez[j][i] * maciez[i][q] / maciez[i][i];
- }
- }
- }
- }
- double metoda_bisekcji(double(*F)(double))
- {
- double a,b;
- double E = 0.000001;
- double i = 0;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- std::cout<< F(a) << " " << F(b) << "nn";
- std::cout << F(a) * F(b) << "nn";
- if (F(a) * F(b) > 0 /*|| (F(a) * F(b) != -INFINITY)*/)
- {
- while (i<100 && F(a) * F(b) > 0)
- {
- b -= b / 100;
- //std::cout << "F(b) = " << F(b) << "n";
- i++;
- }
- if (F(a) * F(b) > 0)
- {
- std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
- return 0.;
- }
- i = 0;
- }
- while(b-a>E && i<60)
- {
- if (F((a + b) / 2) * F(a) < 0)
- b = (a + b) / 2;
- else
- {
- a = (a + b) / 2;
- };
- i++;
- }
- std::cout << "ilosc = " << i <<"nn" ;
- return (a+b)/2;
- }
- double metoda_zgadywana(double(*F)(double))
- {
- double a, b;
- double E = 0.000001;
- double i = 0;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- if (F(a) * F(b) > 0 /*|| (F(a) * F(b) != -INFINITY)*/)
- {
- while (i<100 && F(a) * F(b) > 0)
- {
- b -= b / 100;
- //std::cout << "F(b) = " << F(b) << "n";
- i++;
- }
- if (F(a) * F(b) > 0)
- {
- std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
- return 0.;
- }
- i = 0;
- }
- while (F(a) * F(b) < 0)
- {
- a += (b - a) / 1000;
- }
- return a;
- }
- double zloty_podzial_min(double(*F)(double))
- {
- double a,b;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- double k = (sqrt(5.) - 1.) / 2.;
- double E = 0.00000001;
- int i = 0;
- if (F(a) * F(b) > 0 /*|| (F(a) * F(b) != -INFINITY)*/)
- {
- while (i<100 && F(a) * F(b) > 0)
- {
- b -= b / 100;
- //std::cout << "F(b) = " << F(b) << "n";
- i++;
- }
- if (F(a) * F(b) > 0)
- {
- std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
- return 0.;
- }
- i = 0;
- }
- while (b - a>E && i<60)
- {
- if (F(a + (b - a)*k) * F(a) < 0)
- b = a + (b - a)*k;
- else
- {
- a = a + (b - a)*k;
- };
- i++;
- }
- std::cout << "ilosc = " << i << "nn";
- return (a + b) / 2;
- }
- double zloty_podzial(double(*F)(double))
- {
- double a, b;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- if (F(a) * F(b) > 0)
- {
- std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
- return 0.;
- }
- else
- {
- double k = (sqrt(5.) - 1.) / 2.;
- double xL = b - k * (b - a);
- double xR = a + k * (b - a);
- double E = 0.0000001;
- std::cout << "F(xl) = " << F(xL) << " f(xr) = " << F(xR) << "n";
- std::cout << "xl = " << xL << " xr = " << xR << "n";
- std::cout << "a = " << a << " b = " << b << "n";
- //std::cout << "xl = " << F(xL) << " xr = " << F(xR) << "n";
- while (b - a > E)//dokladność
- {
- std::cout << "F(xl) = " << F(xL) << " f(xr) = " << F(xR) << "n";
- std::cout << "xl = " << xL << " xr = " << xR << "n";
- std::cout << "a = " << a << " b = " << b << "n";
- if (F(xL) < F(xR))
- {
- // wybierz przedział [a, xR]
- b = xR;
- xR = xL;
- xL = b - k * (b - a);
- }
- else
- {
- // wybierz przedział [xL, b]
- a = xL;
- xL = xR;
- xR = a + k * (b - a);
- }
- std::cout << "xl = " << xL << " xr = " << xR << "n";
- }
- return F((a + b) / 2.);
- }
- }
- double metoda_falsi(double(*F)(double))
- {
- double x0, f0, fp;
- double a, b;
- double E = 0.0000001;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- do
- {
- x0 = (a*F(b) - b*F(a)) / (F(b) - F(a));
- f0 = F(x0);
- fp = F(a);
- if (f0 * fp < 0) b = x0;
- else a = x0;
- } while (fabs(f0) > E);
- return x0;
- }
- double metoda_siecznych(double(*F)(double))
- {
- double x0, f0, fp, fk;
- double a, b;
- double E = 0.0000001;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- fp = F(a);
- fk = F(b);
- do
- {
- x0 = (a*fk - b*fp) / (fk - fp);
- f0 = F(x0);
- a = b;
- fp = fk;
- b = x0;
- fk = f0;
- } while (fabs(f0) > E);
- return x0;
- }
- double metoda_newtona_rapshona(double(*F)(double), double(*FP)(double))
- {
- double a, b, x0, x1;
- double E = 0.0000001;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- if (F(a)*F(b) >= 0)
- std::cout << "zly przedzial!!!";
- else
- {
- std::cout << "podaj miejsce startowe: ";
- std::cin >> x0;
- while (fabs(F(x0)) > E )
- {
- x1 = x0 - F(x0) / FP(x0);
- x0 = x1;
- }
- }
- return x0;
- }
- double metoda_newtona_rapshona_2_stopnia(double(*F)(double), double(*FP)(double), double(*FPP)(double))
- {
- double a, b, x0, x1;
- double E = 0.0000001;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- if (F(a)*F(b) >= 0)
- std::cout << "zly przedzial!!!";
- else
- {
- std::cout << "podaj miejsce startowe: ";
- std::cin >> x0;
- while (fabs(F(x0)) > E)
- {
- x1 = x0 - ((FP(x0) - sqrt(pow(FP(x0), 2) - 2 * F(x0)*FPP(x0))) / FPP(x0));
- x0 = x1;
- }
- }
- return x0;
- }
- double metoda_halleya(double(*F)(double), double(*FP)(double), double(*FPP)(double))
- {
- double a, b, x0, x1;
- double E = 0.0000001;
- std::cout << "podaj a :";
- std::cin >> a;
- std::cout << "podaj b :";
- std::cin >> b;
- if (F(a)*F(b) >= 0)
- std::cout << "zly przedzial!!!";
- else
- {
- std::cout << "podaj miejsce startowe: ";
- std::cin >> x0;
- std::cout << x0 << " ";
- while (fabs(F(x0)) > E)
- {
- x1 = x0 - (2 * F(x0)*FP(x0)) / (2 * pow(FP(x0), 2) - F(F(x0))*FPP(x0));
- x0 = x1;
- std::cout << x0 << " ";
- }
- }
- return x0;
- }
- double funkcja(double x)
- {
- if (x == 0)return -1.7e-308;
- else return x - (4 / pow(x, 2));
- }
- double funkcja1(double x)
- {
- if (x == 2)return -1.7e-308;
- else return (pow(x,2)-3)/(x-2);
- }
- double funkcja2(double x)
- {
- return pow(x, 2) - 4 + sin(x);
- }
- double funkcja3(double x)
- {
- return x*x;
- }
- double funkcja4(double x)
- {
- return pow(x,3)+pow(x,2)-3*x-3;
- }
- double funkcja4P(double x)
- {
- return 3*pow(x,2)+2*x-3;
- }
- double funkcja4PP(double x)
- {
- return 6 * x, +2;
- }
- double funkcja5(double x)
- {
- return pow(x, 2) - x + 4;
- }
- double funkcja6(double x)
- {
- return pow(x, 2) + 2*x ;
- }
- double motoda_trapezow(double(*F)(double))
- {
- double a, b, ile,k;
- std::cout << "podaj poczatek przedziału: ";
- std::cin >> a;
- std::cout << "podaj koniec przedziału: ";
- std::cin >> b;
- std::cout << "podaj ile: ";
- std::cin >> ile;
- double suma = 0.;
- k = (b - a) / ile;
- for (int i = 1; i < ile; i++)
- {
- suma += F(a + i * k);
- }
- suma = (suma + (F(a) + F(b)) / 2) * k;
- return suma;
- }
- double metoda_simpsona(double(*F)(double))
- {
- double a, b, ile, k, x,s;
- std::cout << "podaj poczatek przedziału: ";
- std::cin >> a;
- std::cout << "podaj koniec przedziału: ";
- std::cin >> b;
- std::cout << "podaj ile: ";
- std::cin >> ile;
- double suma = 0.;
- s = 0.;
- k = (b - a) / ile;
- for (int i = 1; i <= ile; i++)
- {
- x = a + i*k;
- s += F(x - k / 2);
- if (i < ile) suma += F(x);
- }
- suma = k / 6 * (F(a) + F(b) + 2 * s + 4 * s);
- return suma;
- }
- double metoda_monte_carlo(double(*F)(double))
- {
- srand(time(NULL));
- double a, b, ile, k, x, s;
- std::cout << "podaj poczatek przedziału: ";
- std::cin >> a;
- std::cout << "podaj koniec przedziału: ";
- std::cin >> b;
- std::cout << "podaj ile: ";
- std::cin >> ile;
- double suma = 0.;
- s = 0.;
- k = b - a;
- for (int i = 1; i <= ile; i++)
- {
- suma += F(a + ((double)rand() / (double)(RAND_MAX + 1)*k));
- }
- suma = k * suma / ile;
- return suma;
- }