4回目　課題解答・解説

課題1）

#include <stdio.h>

struct Point2D{
double x;
double y;
};

//ラグランジュ補間の関数
double LagrangeInterp(const Point2D pnt[], const int n, const double x)
/*
* pnt : 点群を表す配列
* n : 点群の個数
* x : 求めたいxの値
* 戻り値 : xに対応する関数値 y
*/
{
　int i,j;
　/* 計算用変数の初期化 */
　double p, s = 0.0;
　for (i = 0; i < n; i++) {
　　p = pnt[i].y; /* 積の計算 */
　　for (j = 0; j < n; j++) {
　　　if (i != j) {
　　　　p *= (x - pnt[j].x) / (pnt[i].x - pnt[j].x);
　　　}
　　}
　　s += p;
　}
　　return s;
}

int main(){
　const int N = 5;
　Point2D points[N] = {
　　{-3.5, -2.2},
　　{-2.2, -0.3},
　　{-1.0, 1.5},
　　{ 0.2, 1.9},
　　{ 2.3, 2.2}
　};
　FILE *fp;

　fp = fopen("Lagrange.csv", "w");
　if(fp == NULL){
　　printf("file cannot be opened.\n");
　　return -1;
　}

　/* xの値を変化させながら LagrangeInterp 関数により y の値を順次計算，出力 */
　for(float x = -4.0; x <= 3.0; x += 0.1){ //-4から3まで計算する
　　fprintf(fp,"%f,%f\n",x,LagrangeInterp(points,N,x));
　}
　fclose(fp);
　printf("出力完了！\n");
　return 0;
}

構造体

関数の設定，引数の表記

constにすることで書き換え不可
構造体初期化，数値代入

ファイルポインタ

課題2）

#include <stdio.h>

struct Point2D
{
　double x;
　double y;
};

const int N = 5; /* データ個数 */

double spline(const Point2D pnt[], const double x)
/*
* スプライン補間関数
* pnt : 既知の点群
* x : xの値
* 戻り値 : xに対応する関数の値
*/
{
　int idx = -1, k;
　double h[N-1], b[N-1], d[N-1], g[N-1], u[N-1], r[N];

　int i;
　for (i=1; i<N-1 && idx<0; i++) {
　　if (x < pnt[i].x )
　　idx = i - 1;
　}

　if (idx < 0)
　　idx = N-2;

　for(i=0; i<N-1; i++) {
　　h[i] = pnt[i+1].x - pnt[i].x;
　}

　for (i=1; i<N-1; i++) {
　　b[i] = 2.0 * (h[i] + h[i-1]);
　　d[i] = 3.0 * ((pnt[i+1].y - pnt[i ].y) / h[i ]
　　- (pnt[i ].y - pnt[i-1].y) / h[i-1]);
　}

　g[1] = h[1] / b[1];

　for (i=2; i<N-2; i++) {
　　g[i] = h[i] / (b[i] - h[i-1] * g[i-1]);
　}

　u[1] = d[1] / b[1];

　for(i=2; i<N-1; i++) {
　　u[i] = (d[i]-h[i-1] * u[i-1]) / (b[i] - h[i-1] * g[i-1]);
　}

　if(idx > 1) {
　　k = idx;
　} else {
　　k = 1;
　}

　r[0] = 0.0;
　r[N-1] = 0.0;
　r[N-2] = u[N-2];

　for (i=N-3; i>=k; i--) {
　　r[i] = u[i] - g[i] * r[i+1];
　}

　double dx = x - pnt[idx].x;
　double q = (pnt[idx+1].y - pnt[idx].y) / h[idx] - h[idx] * (r[idx+1] + 2.0 * r[idx]) / 3.0;
　double s = (r[idx+1] - r[idx]) / (3.0 * h[idx]);

　return pnt[idx].y + dx * (q + dx * (r[idx] + s * dx));
}

int main(){
　Point2D points[N] =
　　{{-3.5, -2.2},
　　{-2.2, -0.3},
　　{-1.0, 1.5},
　　{ 0.2, 1.9},
　　{ 2.3, 2.2}
　};
　FILE *fp;

　fp = fopen("spline.csv", "w");
　if(fp == NULL){
　　printf("file cannot be opened\n");
　　return -1;
　}
　/* spline関数を利用して，あるxの範囲についてyを計算，ファイル出力． */

　for(float x = -4.0; x <= 3.0; x += 0.1){ //-4から3まで計算する
　　fprintf(fp,"%f,%f\n",x,spline(points,x));
　}
　fclose(fp);
　printf("出力完了!\n");
　return 0;
}

構造体

constにより数値を固定

スプライン補間の関数

構造体初期化，数値代入