#include <iostream.h>
#include <iomanip>
#include <math.h>
using namespace std;
int main()
{
cout << 3 + 5 << endl;
cout << 3 - 5 << endl;
cout << 3 * 5 << endl;
cout << pow(3, 5) << endl;
cout << 5 / 3 << endl;
cout << 5.0 / 3 << endl;
cout << 5 / 3.0 << endl;
cout << 5 % 3 << endl;
cout << setw(3) << 3 * 5 << endl;
cout << setw(23) << fixed << setprecision(20) << 5 / 3.0 << endl;
return 0;
}
#include <iostream.h>
int main()
{
int i = 3 * 5;
cout << "3 * 5 = " << i << endl;;
return 0;
}
#include <iostream.h>
int main()
{
for (int i = 1; i < 10; i++)
{
cout << i << ", ";
}
cout << endl;;
return 0;
}
#include <iostream.h>
int main()
{
for (int i = 1; i < 10; i++)
{
if (i % 3 == 0)
{
cout << i << ", ";
}
}
cout << endl;;
return 0;
}
#include <iostream.h>
int main()
{
int sum = 0;
for (int i = 1; i < 100; i++)
{
if (i % 3 == 0)
{
sum += i;
}
}
cout << sum << endl;;
return 0;
}
#include <iostream>
using namespace std;
int sn(int a, int lim);
int main()
{
// 3 の倍数の合計
cout << sn(3, 999) << endl;
return 0;
}
// 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数
int sn(int a, int lim)
{
int n = lim / a; // 項数:n = 上限:lim / 公差:a
int l = n * a; // 末項:l = 項数:n * 公差:a
return (a + l) * n / 2; // 総和:sn = (初項:a + 末項:l) * 項数:n / 2
}
#include <iostream>
using namespace std;
int main()
{
// 10000 までの 自然数の和
// 項目数 n = 10000
int n = 10000;
cout << ( n * (n + 1) / 2 ) << endl;
return 0;
}
#include <iostream>
using namespace std;
int main()
{
// 10000 までの 偶数の和
// 項目数 n = 5000
int n = 10000 / 2;
cout << ( n * (n + 1) ) << endl;
return 0;
}
#include <iostream>
#include <math.h>
using namespace std;
int main()
{
// 10000 までの 奇数の和
// 項目数 n = 5000
int n = 10000 / 2;
cout << int(pow(n, 2)) << endl;
return 0;
}
#include <iostream>
using namespace std;
int main()
{
// 1000 までの 自然数の２乗の和
int n = 1000;
cout << ( n * (n + 1) * (2 * n + 1) / 6 ) << endl;
return 0;
}
#include <iostream>
#include <math.h>
using namespace std;
int main()
{
// 100 までの 自然数の３乗の和
int n = 100;
cout << int(pow(n, 2) * pow((n + 1), 2) / 4 ) << endl;
return 0;
}
#include <iostream>
#include <math.h>
using namespace std;
int main()
{
// 初項 2, 公比 3, 項数 10 の等比数列の和
int n = 10;
int a = 2;
int r = 3;
cout << ( (a * (pow(r, n) - 1)) / (r - 1) ) << endl;
return 0;
}
#include <iostream>
using namespace std;
int main()
{
int a = 5; // 初項 5
int d = 3; // 公差 3
int n = 10; // 項数 10
__int64 p = 1; // 積
for (int i = 1; i <= n; i++)
{
int m = a + (d * (i - 1));
p *= m;
}
cout << p << endl;
return 0;
}
#include <iostream>
using namespace std;
__int64 product(int m, int d, int n);
int main()
{
// 初項 5, 公差 3, 項数 10 の数列の積を表示する
cout << product(5, 3, 10) << endl;
return 0;
}
__int64 product(int m, int d, int n)
{
if (n == 0)
return 1;
else
return m * product(m + d, d, n - 1);
}
#include <iostream>
using namespace std;
// 階乗を求める関数
int Fact(int n)
{
if (n <= 1)
return 1;
else
return n * Fact(n - 1);
}
int main()
{
// 10の階乗
cout << Fact(10) << endl;
cout << (10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1) << endl;
return 0;
}
#include <iostream>
using namespace std;
// 下降階乗冪
int FallingFact(int x, int n)
{
if (n <= 1)
return x;
else
return x * FallingFact(x - 1, n - 1);
}
int main()
{
// 10 から 6 までの 総乗
cout << FallingFact(10, 5) << endl;
cout << (10 * 9 * 8 * 7 * 6) << endl;
return 0;
}
#include <iostream>
using namespace std;
// 上昇階乗冪
int RisingFact(int x, int n)
{
if (n <= 1)
return x;
else
return x * RisingFact(x + 1, n - 1);
}
int main()
{
// 10 から 14 までの 総乗
cout << RisingFact(10, 5) << endl;
cout << (10 * 11 * 12 * 13 * 14) << endl;
return 0;
}
#include <iostream>
using namespace std;
// 階乗
int Fact(int n)
{
if (n <= 1)
return 1;
else
return n * Fact(n - 1);
}
// 下降階乗冪
int FallingFact(int x, int n)
{
if (n <= 1)
return x;
else
return x * FallingFact(x - 1, n - 1);
}
int main()
{
// 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
int n = 10;
int r = 5;
cout << (Fact(n) / Fact(n - r)) << endl;
cout << FallingFact(n, r) << endl;
return 0;
}
#include <iostream>
#include <math.h>
using namespace std;
int main()
{
// 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数)
int n = 10;
int r = 5;
cout << pow(n, r) << endl;
return 0;
}
#include <iostream>
#include <math.h>
using namespace std;
// 組合せ
int Comb(int n, int r)
{
if (r == 0 || r == n)
return 1;
else if (r == 1)
return n;
else
return Comb(n - 1, r - 1) + Comb(n - 1, r);
}
int main()
{
// 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数)
int n = 10;
int r = 5;
cout << Comb(n, r) << endl;
return 0;
}
#include <iostream>
using namespace std;
// 組合せ
int Comb(int n, int r)
{
if (r == 0 || r == n)
return 1;
else if (r == 1)
return n;
else
return Comb(n - 1, r - 1) + Comb(n - 1, r);
}
int main()
{
// 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数)
int n = 10;
int r = 5;
cout << Comb(n + r - 1, r) << endl;
return 0;
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double mySin(double x, int n, bool nega, double numerator, double denominator, double y);
int main()
{
for (int degree = 0; degree <= 360; degree += 15)
{
if (degree % 30 == 0 || degree % 45 == 0)
{
double radian = degree * M_PI / 180.0;
// 自作の正弦関数
double d1 = mySin(radian, 1, false, radian, 1.0, radian);
// 標準の正弦関数
double d2 = sin(radian);
// 標準関数との差異
cout << setw(3) << degree << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
}
return 0;
}
// 自作の正弦関数
double mySin(double x, int n, bool nega, double numerator, double denominator, double y)
{
int m = 2 * n;
denominator = denominator * (m + 1) * m;
numerator = numerator * x * x;
double a = numerator / denominator;
// 十分な精度になったら処理を抜ける
if (a <= 0.00000000001)
return y;
else
return y + mySin(x, ++n, !nega, numerator, denominator, nega ? a : -a);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myCos(double x, int n, bool nega, double numerator, double denominator, double y);
int main()
{
for (int degree = 0; degree <= 360; degree += 15)
{
if (degree % 30 == 0 || degree % 45 == 0)
{
double radian = degree * M_PI / 180.0;
// 自作の余弦関数
double d1 = myCos(radian, 1, false, 1.0, 1.0, 1.0);
// 標準の余弦関数
double d2 = cos(radian);
// 標準関数との差異
cout << setw(3) << degree << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
}
return 0;
}
// 自作の余弦関数
double myCos(double x, int n, bool nega, double numerator, double denominator, double y)
{
int m = 2 * n;
denominator = denominator * m * (m - 1);
numerator = numerator * x * x;
double a = numerator / denominator;
// 十分な精度になったら処理を抜ける
if (a <= 0.00000000001)
return y;
else
return y + myCos(x, ++n, !nega, numerator, denominator, nega ? a : -a);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myTan(double x, double x2, int n, double t);
int main()
{
for (int degree = -90; degree <= 90; degree += 15)
{
if ((degree + 90) % 180 != 0)
{
double radian = degree * M_PI / 180.0;
double x2 = radian * radian;
// 自作の正接関数
double d1 = myTan(radian, x2, 15, 0.0); // 15:必要な精度が得られる十分大きな奇数
// 標準の正接関数
double d2 = tan(radian);
// 標準関数との差異
cout << setw(3) << degree << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
}
return 0;
}
// 自作の正接関数
double myTan(double x, double x2, int n, double t)
{
t = x2 / (n - t);
n -= 2;
if (n <= 1)
return x / (1 - t);
else
return myTan(x, x2, n, t);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myExp(double x, int n, double numerator, double denominator, double y);
int main()
{
for (int i = -10; i <= 10; i++)
{
double x = i / 4.0;
// 標準の指数関数
double d1 = exp(x);
// 自作の指数関数
double d2 = myExp(x, 1, 1.0, 1.0, 1.0);
// 標準関数との差異
cout << setw(5) << fixed << setprecision(2) << x << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
return 0;
}
// 自作の指数関数
double myExp(double x, int n, double numerator, double denominator, double y)
{
denominator = denominator * n;
numerator = numerator * x;
double a = numerator / denominator;
// 十分な精度になったら処理を抜ける
if (fabs(a) <= 0.00000000001)
return y;
else
return y + myExp(x, ++n, numerator, denominator, a);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myExp(double x, double x2, int n, double t);
int main()
{
for (int i = -10; i <= 10; i++)
{
double x = i / 4.0;
// 標準の指数関数
double d1 = exp(x);
// 自作の指数関数
double x2 = x * x;
double d2 = myExp(x, x2, 30, 0.0); // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる
// 標準関数との差異
cout << setw(5) << fixed << setprecision(2) << x << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
return 0;
}
// 自作の指数関数
double myExp(double x, double x2, int n, double t)
{
t = x2 / (n + t);
n -= 4;
if (n < 6)
return 1 + ((2 * x) / (2 - x + t));
else
return myExp(x, x2, n, t);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myLog(double x2, double numerator, double denominator, double y);
int main()
{
for (int i = 1; i <= 20; i++)
{
double x = i / 5.0;
// 標準の対数関数
double d1 = log(x);
// 自作の対数関数
double x2 = (x - 1) / (x + 1);
double d2 = 2 * myLog(x2, x2, 1.0, x2);
// 標準関数との差異
cout << setw(5) << fixed << setprecision(2) << x << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
return 0;
}
// 自作の対数関数
double myLog(double x2, double numerator, double denominator, double y)
{
denominator = denominator + 2;
numerator = numerator * x2 * x2;
double a = numerator / denominator;
// 十分な精度になったら処理を抜ける
if (fabs(a) <= 0.00000000001)
return y;
else
return y + myLog(x2, numerator, denominator, a);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myLog(double x, int n, double t);
int main()
{
for (int i = 1; i <= 20; i++)
{
double x = i / 5.0;
// 標準の対数関数
double d1 = log(x);
// 自作の対数関数
double d2 = myLog(x - 1, 27, 0.0); // 27:必要な精度が得られる十分大きな奇数
// 標準関数との差異
cout << setw(5) << fixed << setprecision(2) << x << ":";
cout << setw(14) << fixed << setprecision(10) << d1 << " - ";
cout << setw(14) << fixed << setprecision(10) << d2 << " = ";
cout << setw(14) << fixed << setprecision(10) << d1 - d2 << endl;
}
return 0;
}
// 自作の対数関数
double myLog(double x, int n, double t)
{
int n2 = n;
double x2 = x;
if (n > 3)
{
if (n % 2 == 0)
n2 = 2;
x2 = x * (n / 2);
}
t = x2 / (n2 + t);
if (n <= 2)
return x / (1 + t);
else
return myLog(x, --n, t);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double mySinh(double x, int n, double numerator, double denominator, double y);
int main()
{
for (int x = -10; x <= 10; x++)
{
// 自作の双曲線正弦関数
double d1 = mySinh(x, 1, x, 1.0, x);
// 標準の双曲線正弦関数
double d2 = sinh(x);
// 標準関数との差異
cout << setw(3) << x << ":";
cout << setw(17) << fixed << setprecision(10) << d1 << " - ";
cout << setw(17) << fixed << setprecision(10) << d2 << " = ";
cout << setw(17) << fixed << setprecision(10) << d1 - d2 << endl;
}
return 0;
}
// 自作の双曲線正弦関数
double mySinh(double x, int n, double numerator, double denominator, double y)
{
int m = 2 * n;
denominator = denominator * (m + 1) * m;
numerator = numerator * x * x;
double a = numerator / denominator;
// 十分な精度になったら処理を抜ける
if (fabs(a) <= 0.00000000001)
return y;
else
return y + mySinh(x, ++n, numerator, denominator, a);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double myCosh(double x, int n, double numerator, double denominator, double y);
int main()
{
for (int x = -10; x <= 10; x++)
{
// 自作の双曲線余弦関数
double d1 = myCosh(x, 1, 1.0, 1.0, 1.0);
// 標準の双曲線余弦関数
double d2 = cosh(x);
// 標準関数との差異
cout << setw(3) << x << ":";
cout << setw(17) << fixed << setprecision(10) << d1 << " - ";
cout << setw(17) << fixed << setprecision(10) << d2 << " = ";
cout << setw(17) << fixed << setprecision(10) << d1 - d2 << endl;
}
return 0;
}
// 自作の双曲線余弦関数
double myCosh(double x, int n, double numerator, double denominator, double y)
{
int m = 2 * n;
denominator = denominator * m * (m - 1);
numerator = numerator * x * x;
double a = numerator / denominator;
// 十分な精度になったら処理を抜ける
if (fabs(a) <= 0.00000000001)
return y;
else
return y + myCosh(x, ++n, numerator, denominator, a);
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double f(double x)
{
return 4 / (1 + x * x);
}
int main()
{
const double a = 0;
const double b = 1;
// 台形則で積分
int n = 2;
for (int j = 1; j <= 10; j++)
{
double h = (b - a) / n;
double s = 0;
double x = a;
for (int i = 1; i <= n - 1; i++)
{
x += h;
s += f(x);
}
s = h * ((f(a) + f(b)) / 2 + s);
n *= 2;
// 結果を π と比較
cout << setw(2) << j << ":";
cout << setw(13) << fixed << setprecision(10) << s << ", ";
cout << setw(13) << fixed << setprecision(10) << s - M_PI << endl;
}
return 0;
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double f(double x)
{
return 4 / (1 + x * x);
}
int main()
{
const double a = 0;
const double b = 1;
// 中点則で積分
int n = 2;
for (int j = 1; j <= 10; j++)
{
double h = (b - a) / n;
double s = 0;
double x = a + (h / 2);
for (int i = 1; i <= n; i++)
{
s += f(x);
x += h;
}
s *= h;
n *= 2;
// 結果を π と比較
cout << setw(2) << j << ":";
cout << setw(13) << fixed << setprecision(10) << s << ", ";
cout << setw(13) << fixed << setprecision(10) << s - M_PI << endl;
}
return 0;
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double f(double x)
{
return 4 / (1 + x * x);
}
int main()
{
const double a = 0;
const double b = 1;
// Simpson則で積分
int n = 2;
for (int j = 1; j <= 5; j++)
{
double h = (b - a) / n;
double s2 = 0;
double s4 = 0;
double x = a + h;
for (int i = 1; i <= n / 2; i++)
{
s4 += f(x);
x += h;
s2 += f(x);
x += h;
}
s2 = (s2 - f(b)) * 2 + f(a) + f(b);
s4 *= 4;
double s = (s2 + s4) * h / 3;
n *= 2;
// 結果を π と比較
cout << setw(2) << j << ":";
cout << setw(13) << fixed << setprecision(10) << s << ", ";
cout << setw(13) << fixed << setprecision(10) << s - M_PI << endl;
}
return 0;
}
#include <iostream>
#include <iomanip>
#include <math>
using namespace std;
double f(double x)
{
return 4 / (1 + x * x);
}
int main()
{
const double a = 0;
const double b = 1;
double t[7][7];
// 台形則で積分
int n = 2;
for (int i = 1; i <= 6; i++)
{
double h = (b - a) / n;
double s = 0;
double x = a;
for (int j = 1; j <= n - 1; j++)
{
x += h;
s += f(x);
}
// 結果を保存
t[i][1] = h * ((f(a) + f(b)) / 2 + s);
n *= 2;
}
// Richardsonの補外法
n = 4;
for (int j = 2; j <= 6; j++)
{
for (int i = j; i <= 6; i++)
{
t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1);
if (i == j)
{
// 結果を π と比較
cout << setw(2) << j << ":";
cout << setw(13) << fixed << setprecision(10) << t[i][j] << ", ";
cout << setw(13) << fixed << setprecision(10) << t[i][j] - M_PI << endl;
}
}
n *= 4;
}
return 0;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// データ点の数
const int N = 7;
// 元の関数
double f(double x);
// Lagrange (ラグランジュ) 補間
double lagrange(double d, double x[], double y[]);
int main()
{
double x[N], y[N];
// 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for (int i = 0; i < N; i++)
{
double d = i * 1.5 - 4.5;
x[i] = d;
y[i] = f(d);
}
// 0.5刻みで 与えられていない値を補間
for (int i = 0; i <= 18; i++)
{
double d = i * 0.5 - 4.5;
double d1 = f(d);
double d2 = lagrange(d, x, y);
// 元の関数と比較
cout << setw(5) << fixed << setprecision(2) << d << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 << '\t';
cout << setw(8) << fixed << setprecision(5) << d2 << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
}
return 0;
}
// 元の関数
double f(double x)
{
return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// Lagrange (ラグランジュ) 補間
double lagrange(double d, double x[], double y[])
{
double sum = 0;
for (int i = 0; i < N; i++)
{
double prod = y[i];
for (int j = 0; j < N; j++)
{
if (j != i)
prod *= (d - x[j]) / (x[i] - x[j]);
}
sum += prod;
}
return sum;
}
#include <iostream>
#include <iomanip>
#include <math.h>
#include <stdio.h>
using namespace std;
// データ点の数
const int N = 7;
// 元の関数
double f(double x);
// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[]);
int main()
{
double x[N], y[N];
// 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for (int i = 0; i < N; i++)
{
double d = i * 1.5 - 4.5;
x[i] = d;
y[i] = f(d);
}
// 0.5刻みで 与えられていない値を補間
for (int i = 0; i <= 18; i++)
{
double d = i * 0.5 - 4.5;
double d1 = f(d);
double d2 = neville(d, x, y);
// 元の関数と比較
cout << setw(5) << fixed << setprecision(2) << d << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 << '\t';
cout << setw(8) << fixed << setprecision(5) << d2 << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
}
return 0;
}
// 元の関数
double f(double x)
{
return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[])
{
double w[N][N];
for (int i = 0; i < N; i++)
w[0][i] = y[i];
for (int j = 1; j < N; j++)
{
for (int i = 0; i < N - j; i++)
w[j][i] = w[j-1][i+1] + (w[j-1][i+1] - w[j-1][i]) * (d - x[i+j]) / (x[i+j] - x[i]);
}
return w[N-1][0];
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// データ点の数
const int N = 7;
// 元の関数
double f(double x);
// Newton (ニュートン) 補間
double newton(double d, double x[], double a[]);
int main()
{
double x[N], y[N];
// 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for (int i = 0; i < N; i++)
{
double d = i * 1.5 - 4.5;
x[i] = d;
y[i] = f(d);
}
// 差分商の表を作る
double d[N][N];
for (int j = 0; j < N; j++)
d[0][j] = y[j];
for (int i = 1; i < N; i++)
{
for (int j = 0; j < N - i; j++)
d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]);
}
// ｎ階差分商
double a[N];
for (int j = 0; j < N; j++)
a[j] = d[j][0];
// 0.5刻みで 与えられていない値を補間
for (int i = 0; i <= 18; i++)
{
double d = i * 0.5 - 4.5;
double d1 = f(d);
double d2 = newton(d, x, a);
// 元の関数と比較
cout << setw(5) << fixed << setprecision(2) << d << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 << '\t';
cout << setw(8) << fixed << setprecision(5) << d2 << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
}
return 0;
}
// 元の関数
double f(double x)
{
return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// Newton (ニュートン) 補間
double newton(double d, double x[], double a[])
{
double sum = a[0];
for (int i = 1; i < N; i++)
{
double prod = a[i];
for (int j = 0; j < i; j++)
prod *= (d - x[j]);
sum += prod;
}
return sum;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// データ点の数
const int N = 7;
const int Nx2 = 14;
// 元の関数
double f(double x);
// 導関数
double fd(double x);
// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[]);
int main()
{
double x[N], y[N], yd[N];
// 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for (int i = 0; i < N; i++)
{
double d = i * 1.5 - 4.5;
x[i] = d;
y[i] = f(d);
yd[i] = fd(d);
}
// 差分商の表を作る
double z[Nx2];
double d[Nx2][Nx2];
for (int i = 0; i < Nx2; i++)
{
int j = i / 2;
z[i] = x[j];
d[0][i] = y[j];
}
for (int i = 1; i < Nx2; i++)
{
for (int j = 0; j < Nx2 - i; j++)
{
if (i == 1 && j % 2 == 0)
d[i][j] = yd[j / 2];
else
d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]);
}
}
// ｎ階差分商
double a[Nx2];
for (int j = 0; j < Nx2; j++)
a[j] = d[j][0];
// 0.5刻みで 与えられていない値を補間
for (int i = 0; i <= 18; i++)
{
double d = i * 0.5 - 4.5;
double d1 = f(d);
double d2 = hermite(d, z, a);
// 元の関数と比較
cout << setw(5) << fixed << setprecision(2) << d << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 << '\t';
cout << setw(8) << fixed << setprecision(5) << d2 << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
}
return 0;
}
// 元の関数
double f(double x)
{
return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// 導関数
double fd(double x)
{
return 1 - pow(x,2) / 2 + pow(x,4) / (4 * 3 * 2);
}
// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[])
{
double sum = a[0];
for (int i = 1; i < Nx2; i++)
{
double prod = a[i];
for (int j = 0; j < i; j++)
prod *= (d - z[j]);
sum += prod;
}
return sum;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// データ点の数
const int N = 7;
// 元の関数
double f(double x);
// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[]);
int main()
{
double x[N], y[N];
// 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for (int i = 0; i < N; i++)
{
double d = i * 1.5 - 4.5;
x[i] = d;
y[i] = f(d);
}
// ３項方程式の係数の表を作る
double a[N], b[N], c[N], d[N];
for (int i = 1; i < N - 1; i++)
{
a[i] = x[i] - x[i-1];
b[i] = 2.0 * (x[i+1] - x[i-1]);
c[i] = x[i+1] - x[i];
d[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / (x[i] - x[i-1]));
}
// ３項方程式を解く (ト−マス法)
double g[N], s[N];
g[1] = b[1];
s[1] = d[1];
for (int i = 2; i < N - 1; i++)
{
g[i] = b[i] - a[i] * c[i-1] / g[i-1];
s[i] = d[i] - a[i] * s[i-1] / g[i-1];
}
double z[N];
z[0] = 0;
z[N-1] = 0;
z[N-2] = s[N-2] / g[N-2];
for (int i = N - 3; i >= 1; i--)
{
z[i] = (s[i] - c[i] * z[i+1]) / g[i];
}
// 0.5刻みで 与えられていない値を補間
for (int i = 0; i <= 18; i++)
{
double d = i * 0.5 - 4.5;
double d1 = f(d);
double d2 = spline(d, x, y, z);
// 元の関数と比較
cout << setw(5) << fixed << setprecision(2) << d << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 << '\t';
cout << setw(8) << fixed << setprecision(5) << d2 << '\t';
cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
}
return 0;
}
// 元の関数
double f(double x)
{
return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[])
{
// 補間関数値がどの区間にあるか
int k = -1;
for (int i = 1; i < N; i++)
{
if (d <= x[i])
{
k = i - 1;
break;
}
}
if (k < 0) k = N - 1;
double d1 = x[k+1] - d;
double d2 = d - x[k];
double d3 = x[k+1] - x[k];
return (z[k] * pow(d1,3) + z[k+1] * pow(d2,3)) / (6.0 * d3)
+ (y[k] / d3 - z[k] * d3 / 6.0) * d1
+ (y[k+1] / d3 - z[k+1] * d3 / 6.0) * d2;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.01;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[2];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[2];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// Euler法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 位置
x += h * vx[0];
y += h * vy[0];
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
// 速度
vx[1] = vx[0] + h * fx(vx[0], vy[0]);
vy[1] = vy[0] + h * fy(vx[0], vy[0]);
vx[0] = vx[1];
vy[0] = vy[1];
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.1;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[2];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[2];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// 中点法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
vx[1] = vx[0] + h / 2 * fx(vx[0], vy[0]);
vy[1] = vy[0] + h / 2 * fy(vx[0], vy[0]);
cout << setw(8) << fixed << setprecision(5) << vx[1] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[1] << "\t";
// 位置
x += h * vx[1];
y += h * vy[1];
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
// 速度
vx[1] = vx[0] + h * fx(vx[0], vy[0]);
vy[1] = vy[0] + h * fy(vx[0], vy[0]);
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.01;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[3];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[3];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x[3];
x[0] = 0.0;
double y[3];
y[0] = 0.0;
// Heun法
for (int i = 1; y[0] >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 位置・速度
x[1] = x[0] + h * vx[0];
y[1] = y[0] + h * vy[0];
vx[1] = vx[0] + h * fx(vx[0], vy[0]);
vy[1] = vy[0] + h * fy(vx[0], vy[0]);
x[2] = x[0] + h * ( vx[0] + vx[1] ) / 2;
y[2] = y[0] + h * ( vy[0] + vy[1] ) / 2;
vx[2] = vx[0] + h * (fx(vx[0], vy[0]) + fx(vx[1], vy[1])) / 2;
vy[2] = vy[0] + h * (fy(vx[0], vy[0]) + fy(vx[1], vy[1])) / 2;
x[0] = x[2];
y[0] = y[2];
vx[0] = vx[2];
vy[0] = vy[2];
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << x[0] << "\t";
cout << setw(8) << fixed << setprecision(5) << y[0] << endl;
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.1;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[3];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[3];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// Heun法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
vx[1] = vx[0] + h * fx(vx[0], vy[0]);
vy[1] = vy[0] + h * fy(vx[0], vy[0]);
vx[2] = (vx[0] + vx[1]) / 2;
vy[2] = (vy[0] + vy[1]) / 2;
cout << setw(8) << fixed << setprecision(5) << vx[2] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[2] << "\t";
// 位置
x += h * vx[2];
y += h * vy[2];
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
// 速度退避
vx[0] = vx[1];
vy[0] = vy[1];
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double dt = 0.1;
// 空気抵抗による水平方向の減速分
double fx(double vx, double v);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[3];
vx[1] = vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[3];
vy[1] = vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// Heun法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * dt;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
vx[2] = vx[1] + dt * fx(vx[1], v);
vy[2] = vy[1] + dt * fy(vy[1], v);
vx[0] = (vx[1] + vx[2]) / 2;
vy[0] = (vy[1] + vy[2]) / 2;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
// 位置
x = x + dt * vx[0];
y = y + dt * vy[0];
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
// 速度退避
vx[1] = vx[2];
vy[1] = vy[2];
v = sqrt(vx[1] * vx[1] + vy[1] * vy[1]);
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double v)
{
return k * v * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v)
{
return g + (k * v * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.01;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[2];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[2];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x[2];
x[0] = 0.0;
double y[2];
y[0] = 0.0;
// 中点法
for (int i = 1; y[0] >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 位置・速度
vx[1] = h * fx(vx[0], vy[0]);
vy[1] = h * fy(vx[0], vy[0]);
double wx = vx[0] + vx[1] / 2;
double wy = vy[0] + vy[1] / 2;
vx[0] = vx[0] + h * fx(wx, wy);
vy[0] = vy[0] + h * fy(wx, wy);
x[0] = x[0] + h * wx;
y[0] = y[0] + h * wy;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << x[0] << "\t";
cout << setw(8) << fixed << setprecision(5) << y[0] << endl;
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.1;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[5];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[5];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// Runge-Kutta法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
vx[1] = h * fx(vx[0], vy[0]);
vy[1] = h * fy(vx[0], vy[0]);
double wx = vx[0] + vx[1] / 2;
double wy = vy[0] + vy[1] / 2;
vx[2] = h * fx(wx, wy);
vy[2] = h * fy(wx, wy);
wx = vx[0] + vx[2] / 2;
wy = vy[0] + vy[2] / 2;
vx[3] = h * fx(wx, wy);
vy[3] = h * fy(wx, wy);
wx = vx[0] + vx[3];
wy = vy[0] + vy[3];
vx[4] = h * fx(wx, wy);
vy[4] = h * fy(wx, wy);
vx[0] += (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6;
vy[0] += (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
// 位置
x += h * vx[0];
y += h * vy[0];
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double dt = 0.1;
// 空気抵抗による水平方向の減速分
double fx(double vx, double v);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[5];
vx[1] = vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[5];
vy[1] = vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// Runge-Kutta法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * dt;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
vx[2] = vx[1] + dt / 2 * fx(vx[1], v) / 2;
vy[2] = vy[1] + dt / 2 * fy(vy[1], v) / 2;
v = sqrt(vx[2] * vx[2] + vy[2] * vy[2]);
vx[3] = vx[2] + dt / 2 * fx(vx[2], v) / 2;
vy[3] = vy[2] + dt / 2 * fy(vy[2], v) / 2;
v = sqrt(vx[3] * vx[3] + vy[3] * vy[3]);
vx[4] = vx[3] + dt * fx(vx[3], v);
vy[4] = vy[3] + dt * fy(vy[3], v);
vx[0] = (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6;
vy[0] = (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
// 位置
x = x + dt * vx[0];
y = y + dt * vy[0];
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
// 速度
v = sqrt(vx[1] * vx[1] + vy[1] * vy[1]);
vx[1] = vx[1] + dt * fx(vx[1], v);
vy[1] = vy[1] + dt * fy(vy[1], v);
v = sqrt(vx[1] * vx[1] + vy[1] * vy[1]);
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double v)
{
return k * v * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v)
{
return g + (k * v * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.01;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[5];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[5];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x[5];
x[0] = 0.0;
double y[5];
y[0] = 0.0;
// Runge-Kutta法
for (int i = 1; y[0] >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 位置・速度
x[1] = h * vx[0];
y[1] = h * vy[0];
vx[1] = h * fx(vx[0], vy[0]);
vy[1] = h * fy(vx[0], vy[0]);
double wx = vx[0] + vx[1] / 2;
double wy = vy[0] + vy[1] / 2;
x[2] = h * wx;
y[2] = h * wy;
vx[2] = h * fx(wx, wy);
vy[2] = h * fy(wx, wy);
wx = vx[0] + vx[2] / 2;
wy = vy[0] + vy[2] / 2;
x[3] = h * wx;
y[3] = h * wy;
vx[3] = h * fx(wx, wy);
vy[3] = h * fy(wx, wy);
wx = vx[0] + vx[3];
wy = vy[0] + vy[3];
x[4] = h * wx;
y[4] = h * wy;
vx[4] = h * fx(wx, wy);
vy[4] = h * fy(wx, wy);
x[0] += ( x[1] + x[2] * 2 + x[3] * 2 + x[4]) / 6;
y[0] += ( y[1] + y[2] * 2 + y[3] * 2 + y[4]) / 6;
vx[0] += (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6;
vy[0] += (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << x[0] << "\t";
cout << setw(8) << fixed << setprecision(5) << y[0] << endl;
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double dt = 0.1;
// 空気抵抗による水平方向の減速分
double fx(double vx, double v);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v[5];
v[1] = v[0] = 250 * 1000 / 3600;
// 水平方向の速度
double vx[5];
vx[1] = vx[0] = v[0] * cos(radian);
// 鉛直方向の速度
double vy[5];
vy[1] = vy[0] = v[0] * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// Runge-Kutta-Gill法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * dt;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
vx[2] = vx[1] + dt / 2 * fx(vx[1], v[1]) / 2;
vy[2] = vy[1] + dt / 2 * fy(vy[1], v[1]) / 2;
v[2] = sqrt(vx[2] * vx[2] + vy[2] * vy[2]);
vx[3] = vx[2] + dt / 2 * (fx(vx[1], v[1]) * ((sqrt(2.0) - 1) / 2) + fx(vx[2], v[2]) * (1 - 1 / sqrt(2.0)));
vy[3] = vy[2] + dt / 2 * (fy(vy[1], v[1]) * ((sqrt(2.0) - 1) / 2) + fy(vy[2], v[2]) * (1 - 1 / sqrt(2.0)));
v[3] = sqrt(vx[3] * vx[3] + vy[3] * vy[3]);
vx[4] = vx[3] + dt * (fx(vx[2], v[2]) * (1 / sqrt(2.0)) + fx(vx[3], v[3]) * (1 + 1 / sqrt(2.0)));
vy[4] = vy[3] + dt * (fy(vy[2], v[2]) * (1 / sqrt(2.0)) + fy(vy[3], v[3]) * (1 + 1 / sqrt(2.0)));
vx[0] = (vx[1] + vx[2] * (2 - sqrt(2.0)) + vx[3] * (2 + sqrt(2.0) ) + vx[4]) / 6;
vy[0] = (vy[1] + vy[2] * (2 - sqrt(2.0)) + vy[3] * (2 + sqrt(2.0) ) + vy[4]) / 6;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
// 位置
x = x + dt * vx[0];
y = y + dt * vy[0];
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
// 速度
v[1] = sqrt(vx[1] * vx[1] + vy[1] * vy[1]);
vx[1] = vx[1] + dt * fx(vx[1], v[1]);
vy[1] = vy[1] + dt * fy(vy[1], v[1]);
v[1] = sqrt(vx[1] * vx[1] + vy[1] * vy[1]);
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double v)
{
return k * v * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v)
{
return g + (k * v * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.01;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[5];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[5];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x[5];
x[0] = 0.0;
double y[5];
y[0] = 0.0;
// Runge-Kutta-Gill法
for (int i = 1; y[0] >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 速度
x[1] = h * vx[0];
y[1] = h * vy[0];
vx[1] = h * fx(vx[0], vy[0]);
vy[1] = h * fy(vx[0], vy[0]);
double wx = vx[0] + vx[1] / 2;
double wy = vy[0] + vy[1] / 2;
x[2] = h * wx;
y[2] = h * wy;
vx[2] = h * fx(wx, wy);
vy[2] = h * fy(wx, wy);
wx = vx[0] + vx[1] * ((sqrt(2.0) - 1) / 2) + vx[2] * (1 - 1 / sqrt(2.0));
wy = vy[0] + vy[1] * ((sqrt(2.0) - 1) / 2) + vy[2] * (1 - 1 / sqrt(2.0));
x[3] = h * wx;
y[3] = h * wy;
vx[3] = h * fx(wx, wy);
vy[3] = h * fy(wx, wy);
wx = vx[0] - vx[2] / sqrt(2.0) + vx[3] * (1 + 1 / sqrt(2.0));
wy = vy[0] - vy[2] / sqrt(2.0) + vy[3] * (1 + 1 / sqrt(2.0));
x[4] = h * wx;
y[4] = h * wy;
vx[4] = h * fx(wx, wy);
vy[4] = h * fy(wx, wy);
x[0] += ( x[1] + x[2] * (2 - sqrt(2.0)) + x[3] * (2 + sqrt(2.0)) + x[4]) / 6;
y[0] += ( y[1] + y[2] * (2 - sqrt(2.0)) + y[3] * (2 + sqrt(2.0)) + y[4]) / 6;
vx[0] += (vx[1] + vx[2] * (2 - sqrt(2.0)) + vx[3] * (2 + sqrt(2.0)) + vx[4]) / 6;
vy[0] += (vy[1] + vy[2] * (2 - sqrt(2.0)) + vy[3] * (2 + sqrt(2.0)) + vy[4]) / 6;
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << x[0] << "\t";
cout << setw(8) << fixed << setprecision(5) << y[0] << endl;
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double h = 0.01;
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy);
int main()
{
// 角度
double degree = 45;
double radian = degree * M_PI / 180.0;
// 初速 250 km/h -> 秒速に変換
double v = 250 * 1000 / 3600;
// 水平方向の速度
double vx[3];
vx[0] = v * cos(radian);
// 鉛直方向の速度
double vy[3];
vy[0] = v * sin(radian);
// 経過秒数
double t = 0.0;
// 位置
double x = 0.0;
double y = 0.0;
// 後退 Euler法
for (int i = 1; y >= 0.0; i++)
{
// 経過秒数
t = i * h;
cout << setw(4) << fixed << setprecision(2) << t << "\t";
// 位置・速度
vx[1] = vx[0] + h * fx(vx[0], vy[0]);
vy[1] = vy[0] + h * fy(vx[0], vy[0]);
for (int j = 0; j <= 10; j++)
{
vx[2] = vx[0] + h * fx(vx[1], vy[1]);
vy[2] = vy[0] + h * fy(vx[1], vy[1]);
if ((fabs(vx[1] - vx[2]) < 0.00001) && (fabs(vy[1] - vy[2]) < 0.00001)) break;
vx[1] = vx[2];
vy[1] = vy[2];
}
vx[0] = vx[1];
vy[0] = vy[1];
x += h * vx[0];
y += h * vy[0];
cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t";
cout << setw(9) << fixed << setprecision(5) << x << "\t";
cout << setw(8) << fixed << setprecision(5) << y << endl;
}
return 0;
}
// 空気抵抗による水平方向の減速分
double fx(double vx, double vy)
{
return k * sqrt(vx * vx + vy * vy) * vx;
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vx, double vy)
{
return g + (k * sqrt(vx * vx + vy * vy) * vy);
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double f(double x)
{
return x * x - 2;
}
double bisection(double a, double b)
{
double c;
while (true)
{
// 区間 (a, b) の中点 c = (a + b) / 2
c = (a + b) / 2;
cout << setw(12) << fixed << setprecision(10) << c << "\t";
cout << setw(13) << fixed << setprecision(10) << c - sqrt(2) << endl;
double fc = f(c);
if (fabs(fc) < 0.0000000001) break;
if (fc < 0)
{
// f(c) < 0 であれば, 解は区間 (c, b) の中に存在
a = c;
}
else
{
// f(c) > 0 であれば, 解は区間 (a, c) の中に存在
b = c;
}
}
return c;
}
int main()
{
double a = 1;
double b = 2;
cout << setw(12) << fixed << setprecision(10) << bisection(a, b) << endl;
return 0;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double f(double x)
{
return x * x - 2;
}
double falseposition(double a, double b)
{
double c;
while (true)
{
// 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点
c = (a * f(b) - b * f(a)) / (f(b) - f(a));
cout << setw(12) << fixed << setprecision(10) << c << "\t";
cout << setw(13) << fixed << setprecision(10) << c - sqrt(2) << endl;
double fc = f(c);
if (fabs(fc) < 0.0000000001) break;
if (fc < 0)
{
// f(c) < 0 であれば, 解は区間 (c, b) の中に存在
a = c;
}
else
{
// f(c) > 0 であれば, 解は区間 (a, c) の中に存在
b = c;
}
}
return c;
}
int main()
{
double a = 1;
double b = 2;
cout << setw(12) << fixed << setprecision(10) << falseposition(a, b) << endl;
return 0;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double g(double x)
{
return (x / 2) + (1 / x);
}
double iterative(double x0)
{
double x1;
while (true)
{
x1 = g(x0);
cout << setw(12) << fixed << setprecision(10) << x1 << "\t";
cout << setw(13) << fixed << setprecision(10) << x1 - sqrt(2) << endl;
if (fabs(x1 - x0) < 0.0000000001) break;
x0 = x1;
}
return x1;
}
int main()
{
double x = 1.0;
cout << setw(12) << fixed << setprecision(10) << iterative(x) << endl;
return 0;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double f0(double x)
{
return x * x - 2;
}
double f1(double x)
{
return 2 * x;
}
double newton(double x0)
{
double x1;
while (true)
{
x1 = x0 - (f0(x0) / f1(x0));
cout << setw(12) << fixed << setprecision(10) << x1 << "\t";
cout << setw(13) << fixed << setprecision(10) << x1 - sqrt(2) << endl;
if (fabs(x1 - x0) < 0.0000000001) break;
x0 = x1;
}
return x1;
}
int main()
{
double x = 2.0;
cout << setw(12) << fixed << setprecision(10) << newton(x) << endl;
return 0;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double f0(double x)
{
return x * x - 2;
}
double f1(double x)
{
return 2 * x;
}
double f2(double x)
{
return 2;
}
double bailey(double x0)
{
double x1;
while (true)
{
x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
cout << setw(12) << fixed << setprecision(10) << x1 << "\t";
cout << setw(13) << fixed << setprecision(10) << x1 - sqrt(2) << endl;
if (fabs(x1 - x0) < 0.0000000001) break;
x0 = x1;
}
return x1;
}
int main()
{
double x = 2.0;
cout << setw(12) << fixed << setprecision(10) << bailey(x) << endl;
return 0;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double f(double x)
{
return x * x - 2;
}
double secant(double x0, double x1)
{
double x2;
while (true)
{
x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0));
cout << setw(12) << fixed << setprecision(10) << x2 << "\t";
cout << setw(13) << fixed << setprecision(10) << x2 - sqrt(2) << endl;
if (fabs(x2 - x1) < 0.0000000001) break;
x0 = x1;
x1 = x2;
}
return x2;
}
int main()
{
double x0 = 1;
double x1 = 2;
cout << setw(12) << fixed << setprecision(10) << secant(x0, x1) << endl;
return 0;
}
#include <iostream>
using namespace std;
void multiply(double a[3][5], double b[5][4], double c[3][4]);
int main()
{
cout << 'A' << endl;
double a[4][4] = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}};
for (int i = 0;i<=3;i++)
{
for (int j = 0;j<=3;j++)
cout << a[i][j] << '\t';
cout << endl;
}
cout << 'B' << endl;
double b[4] = {20,16,8,17};
for (int i = 0;i<=3;i++)
cout << b[i] << '\t';
}
cout << endl;
cout << 'C' << endl;
double c[4]={0,0,0,0};
double r[4]={0,0,0,0};
calc(a,b,c,r);
for (int i = 0;i<=3;i++)
{
cout << c[i] << '\t';
}
cout << endl;
return 0;
}
void calc(double a[4][4], double b[4], double c[4], double r[4])
{
double eps = 0.0000000001
int lst = 1000
double ErrNum;
double xx[4];
for (int l = 0;l<=lst;l++)
{
ErrNum = 0;
for (int i = 0;i<=3;i++)
{
xx[i] = b[i];
for (int j = 0;j<=(i - 1);j++)
xx[i] = xx[i] - a[i][j] * c[j];
for (int j = i + 1;j<=3;j++)
xx[i] = xx[i] - a[i][j] * c[j];
xx[i] = xx[i] / a[i][i];
errnum += fabs(xx[i] - x[i]);
}
for (int i = 0;i<=3;i++)
c[i] = xx[i];
if (errnum <= eps)
{
for (int i = 0;i<=3;i++)
{
r[i] = 0;
for (int j = 0;j<=3;j++)
r[i] += a[i][j] * c[j];
r[i] = r[i] - b[i];
}
return;
}
}
puts("収束しない");
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// ヤコビの反復法
void jacobi(double a[N][N], double b[N], double c[N]);
int main()
{
double a[N][N] = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}};
double b[N] = {20,16,8,17};
double c[N] = {0,0,0,0};
// ヤコビの反復法
jacobi(a,b,c);
cout << "解" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << c[i] << "\t";
cout << endl;
return 0;
}
// ヤコビの反復法
void jacobi(double a[N][N], double b[N], double x0[N])
{
while (true)
{
double x1[N];
bool finish = true;
for (int i = 0; i < N; i++)
{
x1[i] = 0;
for (int j = 0; j < N; j++)
if (j != i)
x1[i] += a[i][j] * x0[j];
x1[i] = (b[i] - x1[i]) / a[i][i];
if (fabs(x1[i] - x0[i]) > 0.0000000001) finish = false;
}
for (int i = 0; i < N; i++)
{
x0[i] = x1[i];
cout << setw(14) << fixed << setprecision(10) << x0[i] << "\t";
}
cout << endl;
if (finish) return;
}
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// ガウス・ザイデル法
void gauss(double a[N][N], double b[N], double c[N]);
int main()
{
double a[N][N] = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}};
double b[N] = {20,16,8,17};
double c[N] = {0,0,0,0};
// ガウス・ザイデル法
gauss(a,b,c);
cout << "解" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << c[i] << "\t";
cout << endl;
return 0;
}
// ガウス・ザイデル法
void gauss(double a[N][N], double b[N], double x0[N])
{
while (true)
{
double x1;
bool finish = true;
for (int i = 0; i < N; i++)
{
x1 = 0;
for (int j = 0; j < N; j++)
if (j != i)
x1 += a[i][j] * x0[j];
x1 = (b[i] - x1) / a[i][i];
if (fabs(x1 - x0[i]) > 0.0000000001) finish = false;
x0[i] = x1;
}
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << x0[i] << "\t";
cout << endl;
if (finish) return;
}
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// ピボット選択
void pivoting(double a[N][N], double b[N]);
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// 後退代入
void backward_substitution(double a[N][N], double b[N]);
// 状態を確認
void disp_progress(double a[N][N], double b[N]);
// ガウスの消去法
int main()
{
double a[N][N] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
double b[N] = {8,17,20,16};
// ピボット選択
pivoting(a,b);
cout << "ピボット選択後" << endl;
disp_progress(a,b);
// 前進消去
forward_elimination(a,b);
cout << "前進消去後" << endl;
disp_progress(a,b);
// 後退代入
backward_substitution(a,b);
cout << "解" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << b[i] << "\t";
cout << endl;
return 0;
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
for (int pivot = 0; pivot < N - 1; pivot++)
{
for (int row = pivot + 1; row < N; row++)
{
double s = a[row][pivot] / a[pivot][pivot];
for (int col = pivot; col < N; col++)
a[row][col] -= a[pivot][col] * s;
b[row] -= b[pivot] * s;
}
}
}
// 後退代入
void backward_substitution(double a[N][N], double b[N])
{
for (int row = N - 1; row >= 0; row--)
{
for (int col = N - 1; col > row; col--)
b[row] -= a[row][col] * b[col];
b[row] /= a[row][row];
}
}
// ピボット選択
void pivoting(double a[N][N], double b[N])
{
for(int pivot = 0; pivot < N; pivot++)
{
// 各列で 一番値が大きい行を 探す
int max_row = pivot;
double max_val = 0;
for (int row = pivot; row < N; row++)
{
if (fabs(a[row][pivot]) > max_val)
{
// 一番値が大きい行
max_val = fabs(a[row][pivot]);
max_row = row;
}
}
// 一番値が大きい行と入れ替え
if (max_row != pivot)
{
double tmp;
for (int col = 0; col < N; col++)
{
tmp = a[max_row][col];
a[max_row][col] = a[pivot][col];
a[pivot][col] = tmp;
}
tmp = b[max_row];
b[max_row] = b[pivot];
b[pivot] = tmp;
}
}
}
// 状態を確認
void disp_progress(double a[N][N], double b[N])
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
cout << setw(14) << fixed << setprecision(10) << b[i] << endl;
}
cout << endl;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// ピボット選択
void pivoting(double a[N][N], double b[N]);
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// 後退代入
void backward_substitution(double a[N][N], double b[N]);
// 状態を確認
void disp_progress(double a[N][N], double b[N]);
// ガウス・ジョルダン法
int main()
{
double a[N][N] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
double b[N] = {8,17,20,16};
// ピボット選択
pivoting(a,b);
cout << "ピボット選択後" << endl;
disp_progress(a,b);
// 前進消去
forward_elimination(a,b);
cout << "前進消去後" << endl;
disp_progress(a,b);
// 後退代入
backward_substitution(a,b);
cout << "解" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << b[i] << "\t";
cout << endl;
return 0;
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
// 対角線上の係数を 1 にする
for (int pivot = 0; pivot < N; pivot++)
{
double s = a[pivot][pivot];
for (int col = 0; col < N; col++)
a[pivot][col] /= s;
b[pivot] /= s;
}
cout << "対角線上の係数を 1 にする" << endl;
disp_progress(a,b);
// 対角行列にする
for (int pivot = 0; pivot < N; pivot++)
{
for (int row = 0; row < N; row++)
{
if (row == pivot) continue;
double s = a[row][pivot] / a[pivot][pivot];
for (int col = pivot; col < N; col++)
a[row][col] -= a[pivot][col] * s;
b[row] -= b[pivot] * s;
}
}
}
// 後退代入
void backward_substitution(double a[N][N], double b[N])
{
for (int pivot = 0; pivot < N; pivot++)
b[pivot] /= a[pivot][pivot];
}
// ピボット選択
void pivoting(double a[N][N], double b[N])
{
for(int pivot = 0; pivot < N; pivot++)
{
// 各列で 一番値が大きい行を 探す
int max_row = pivot;
double max_val = 0;
for (int row = pivot; row < N; row++)
{
if (fabs(a[row][pivot]) > max_val)
{
// 一番値が大きい行
max_val = fabs(a[row][pivot]);
max_row = row;
}
}
// 一番値が大きい行と入れ替え
if (max_row != pivot)
{
double tmp;
for (int col = 0; col < N; col++)
{
tmp = a[max_row][col];
a[max_row][col] = a[pivot][col];
a[pivot][col] = tmp;
}
tmp = b[max_row];
b[max_row] = b[pivot];
b[pivot] = tmp;
}
}
}
// 状態を確認
void disp_progress(double a[N][N], double b[N])
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
cout << setw(14) << fixed << setprecision(10) << b[i] << endl;
}
cout << endl;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// ピボット選択
void pivoting(double a[N][N], double b[N]);
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// LU分解
void decomp(double a[N][N], double b[N], double x[N]);
// 状態を確認
void disp_progress(double a[N][N], double b[N]);
// LU分解
int main()
{
double a[N][N] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
double b[N] = {8,17,20,16};
// ピボット選択
pivoting(a,b);
cout << "ピボット選択後" << endl;
disp_progress(a,b);
// LU分解
double x[N] = {0,0,0,0};
decomp(a,b,x);
cout << "X" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << x[i] << "\t";
cout << endl;
return 0;
}
// LU分解
void decomp(double a[N][N], double b[N], double x[N])
{
// 前進消去
forward_elimination(a,b);
// 下三角行列を確認
cout << "L" << endl;
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
if (i > j)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
else if (i == j)
cout << setw(14) << fixed << setprecision(10) << 1.0 << "\t";
else
cout << setw(14) << fixed << setprecision(10) << 0.0 << "\t";
cout << endl;
}
cout << endl;
// 上三角行列を確認
cout << "U" << endl;
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
if (i <= j)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
else
cout << setw(14) << fixed << setprecision(10) << 0.0 << "\t";
cout << endl;
}
cout << endl;
// Ly=b から y を求める (前進代入)
double y[N] = {0,0,0,0};
for (int row = 0; row < N; row++)
{
for (int col = 0; col < row; col++)
b[row] -= a[row][col] * y[col];
y[row] = b[row];
}
// y の 値 を確認
cout << "Y" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << y[i] << "\t";
cout << endl << endl;
// Ux=y から x を求める (後退代入)
for (int row = N - 1; row >= 0; row--)
{
for (int col = N - 1; col > row; col--)
y[row] -= a[row][col] * x[col];
x[row] = y[row] / a[row][row];
}
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
for (int pivot = 0; pivot < N - 1; pivot++)
{
for (int row = pivot + 1; row < N; row++)
{
double s = a[row][pivot] / a[pivot][pivot];
for (int col = pivot; col < N; col++)
a[row][col] -= a[pivot][col] * s; // これが 上三角行列
a[row][pivot] = s; // これが 下三角行列
// b[row] -= b[pivot] * s; // この値は変更しない
}
}
}
// ピボット選択
void pivoting(double a[N][N], double b[N])
{
for(int pivot = 0; pivot < N; pivot++)
{
// 各列で 一番値が大きい行を 探す
int max_row = pivot;
double max_val = 0;
for (int row = pivot; row < N; row++)
{
if (fabs(a[row][pivot]) > max_val)
{
// 一番値が大きい行
max_val = fabs(a[row][pivot]);
max_row = row;
}
}
// 一番値が大きい行と入れ替え
if (max_row != pivot)
{
double tmp;
for (int col = 0; col < N; col++)
{
tmp = a[max_row][col];
a[max_row][col] = a[pivot][col];
a[pivot][col] = tmp;
}
tmp = b[max_row];
b[max_row] = b[pivot];
b[pivot] = tmp;
}
}
}
// 状態を確認
void disp_progress(double a[N][N], double b[N])
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
cout << setw(14) << fixed << setprecision(10) << b[i] << endl;
}
cout << endl;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// LL^T分解
void decomp(double a[N][N], double b[N], double x[N]);
// 状態を確認
void disp_progress(double a[N][N]);
// コレスキー法
int main()
{
double a[N][N] = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}};
double b[N] = {34,68,96,125};
cout << "A" << endl;
disp_progress(a);
cout << endl;
// LL^T分解
double x[N] = {0,0,0,0};
decomp(a,b,x);
cout << "X" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << x[i] << "\t";
cout << endl;
return 0;
}
// LL^T分解
void decomp(double a[N][N], double b[N], double x[N])
{
// 前進消去
forward_elimination(a,b);
cout << "L と L^T" << endl;
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
if (j <= i)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
else
cout << setw(14) << fixed << setprecision(10) << a[j][i] << "\t";
cout << endl;
}
cout << endl << endl;
// Ly=b から y を求める (前進代入)
double y[N] = {0,0,0,0};
for (int row = 0; row < N; row++)
{
for (int col = 0; col < row; col++)
b[row] -= a[row][col] * y[col];
y[row] = b[row] / a[row][row];
}
// y の 値 を確認
cout << "Y" << endl;
for (int i = 0; i < N;i++)
cout << setw(14) << fixed << setprecision(10) << y[i] << "\t";
cout << endl << endl;
// Ux=y から x を求める (後退代入)
for (int row = N - 1; row >= 0; row--)
{
for (int col = N - 1; col > row; col--)
y[row] -= a[col][row] * x[col];
x[row] = y[row] / a[row][row];
}
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
for (int pivot = 0; pivot < N; pivot++)
{
double s = 0;
for (int col = 0; col < pivot; col++)
s += a[pivot][col] * a[pivot][col];
// ここで根号の中が負の値になると計算できない！
a[pivot][pivot] = sqrt(a[pivot][pivot] - s);
for (int row = pivot + 1; row < N; row++)
{
s = 0;
for (int col = 0; col < pivot; col++)
s += a[row][col] * a[pivot][col];
a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot];
}
}
}
// 状態を確認
void disp_progress(double a[N][N])
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
cout << endl;
}
cout << endl;
}
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
const int N = 4;
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// LDL^T分解
void decomp(double a[N][N], double b[N], double x[N]);
// 状態を確認
void disp_progress(double a[N][N]);
// 修正コレスキー法
int main()
{
double a[N][N] = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}};
double b[N] = {34,68,96,125};
cout << "A" << endl;
disp_progress(a);
cout << endl;
// LDL^T分解
double x[N] = {0,0,0,0};
decomp(a,b,x);
cout << "X" << endl;
for (int i = 0; i < N; i++)
cout << setw(14) << fixed << setprecision(10) << x[i] << "\t";
cout << endl;
return 0;
}
// LDL^T分解
void decomp(double a[N][N], double b[N], double x[N])
{
// 前進消去
forward_elimination(a,b);
cout << "L と D" << endl;
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
if (j <= i)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
else
cout << setw(14) << fixed << setprecision(10) << a[j][i] << "\t";
cout << endl;
}
cout << endl << endl;
// Ly=b から y を求める (前進代入)
double y[N] = {0,0,0,0};
for (int row = 0; row < N; row++)
{
for (int col = 0; col < row; col++)
b[row] -= a[row][col] * y[col];
y[row] = b[row];
}
// y の 値 を確認
cout << "Y" << endl;
for (int i = 0; i < N;i++)
cout << setw(14) << fixed << setprecision(10) << y[i] << "\t";
cout << endl << endl;
// DL^Tx=y から x を求める (後退代入)
for (int row = N - 1; row >= 0; row--)
{
for (int col = N - 1; col > row; col--)
y[row] -= a[col][row] * a[row][row] * x[col];
x[row] = y[row] / a[row][row];
}
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
for (int pivot = 0; pivot < N; pivot++)
{
double s;
// pivot < k の場合
for (int col = 0; col < pivot; col++)
{
s = a[pivot][col];
for (int k = 0; k < col; k++)
s -= a[pivot][k] * a[col][k] * a[k][k];
a[pivot][col] = s / a[col][col];
}
// pivot == k の場合
s = a[pivot][pivot];
for (int k = 0; k < pivot; k++)
s -= a[pivot][k] * a[pivot][k] * a[k][k];
a[pivot][pivot] = s;
}
}
// 状態を確認
void disp_progress(double a[N][N])
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
cout << setw(14) << fixed << setprecision(10) << a[i][j] << "\t";
cout << endl;
}
cout << endl;
}
|