プログラマ専用SNS ミクプラ

数値計算法の課題として出された問題なのですが、よくわかりません。
8個のデータ点(xi、yi)が与えられたとき、これを近似する3次関数「y= p0+p1x+ p2x^2+p3x^3」を与える。

・正規方程式を解くプログラムを作成し、3次関数の式を求めよ。(少数以下第3位を四捨五入して小数点以下第2までの値で示すこと)
・最小二乗近似のプログラムを作成し、-5から+7まで0.2刻みのデータ点について近似値を求め、グラフを作成せよ。

という問題です。
他に、正規方程式を書き下して係数を求める問題もあり、それは教科書をみればできたのですが、それ以降の問題がどうすればいいのかわかりません。よろしくお願い致します。

データ点は
i : 1 2 3 4 5 6 7 8
xi : -4.0, -2.0, -1.0, 0.0, 1.0, 3.0, 4.0, 6.0
yi : -35.1, 15.1, 15.9, 8.9, 0.1, 0.1, 21.1, 135.0

一応係数を求めるプログラムです。

コード:

#include <stdio.h>
#include <math.h>

#define N 8
#define M 3
#define EPS .0001

double a[M+1][M+2];

int Jordan(void);

int main( int argc, char **argv ){
 double x[N]={-4.0,-2.0,-1.0,0.0,1.0,3.0,4.0,6.0};
 double y[N]={-35.1,15.1,15.9,8.9,0.1,21.1,135.0};
 
 int i,j,k;
 
 for(i=0;i<=M;i++)
 for(j=0;j<=M+1;j++)
 for(k=0;k<=N;k++)
 a[i][j]=0.0;
 
 for(i=0;i<=M;i++)
 for(j=0;j<=M+1;j++)
 for(k=0;k<=N;k++)
 a[j][i] +=pow(x[k],(double)(i+j));
 
 for(j=0;j<=M;j++)
 for(k=0;k<=N;k++)
 a[j][M+1]+=y[k]*pow(x[k],j);
 
 if(Jordan()==1) return 1;
 
 for(i=0;i<=M;i++)
 printf("p%2d = %7.3lf\n",i,a[i][M+1]);
 
 return 0;
}

int Jordan(void)
{
 double pivot,delta;
 int i,j,k;
 for(i=0;i<=M;i++){
 pivot=a[i][i];
 if(fabs( pivot ) < EPS ){
 printf( "ピボットが許容範囲誤差以下\n" );
 return 1;
 }
 for( j=0; j<=M+1; j++ )
 a[i][j] /= pivot;
 
 for( k=0; k<=M; k++ ){
 delta = a[k][i];
 for(j=0; j<=M+1; j++)
 if( k != i )
 a[k][j] -= delta * a[i][j];
 }
 }
return 0;

}

A.M さんが書きました: ↑
6年前
一応係数を求めるプログラムです。

そのプログラムには 3個所間違いがあります。
・double y[N] の初期値が 7個しかありません。
・2個所の for(k=0;k<=N;k++) の k<=N は k<N の間違い。

修正してみました。

コード:

#include <stdio.h>
#include <math.h>

#define N   8
#define M   3
#define EPS 0.0001

double a[M+1][M+2];

int Jordan(void);

int main(void)
{
	double x[N] = {  -4.0, -2.0, -1.0, 0.0, 1.0, 3.0,  4.0,   6.0 };
	double y[N] = { -35.1, 15.1, 15.9, 8.9, 0.1, 0.1, 21.1, 135.0 };
	int i, j, k;

	for (i = 0; i <= M; i++)
		for (j = 0; j <= M + 1; j++)
			for (k = 0; k <= N; k++)
				a[i][j] = 0.0;

	for (i = 0; i <= M; i++)
		for (j = 0; j <= M + 1; j++)
			for (k = 0; k < N; k++)
				a[j][i] += pow(x[k], i + j);

	for (j = 0; j <= M; j++)
		for (k = 0; k < N; k++)
			a[j][M + 1] += y[k] * pow(x[k], j);

	if (Jordan() == 1) return 1;

	for (i = 0; i <= M; i++)
		printf("p%d = %7.3lf\n", i, a[i][M + 1]);

	return 0;
}

int Jordan(void)
{
	double pivot, delta;
	int i, j, k;
	for (i = 0; i <= M; i++) {
		pivot = a[i][i];
		if (fabs(pivot) < EPS) {
			printf("ピボットが許容範囲誤差以下\n");
			return 1;
		}
		for (j = 0; j <= M + 1; j++)
			a[i][j] /= pivot;

		for (k = 0; k <= M; k++) {
			delta = a[k][i];
			for (j = 0; j <= M + 1; j++)
				if (k != i)
					a[k][j] -= delta * a[i][j];
		}
	}
	return 0;
}

実行結果

コード:

p0 =   9.011
p1 =  -8.966
p2 =  -1.000
p3 =   0.999

DxLib を使って図を書くプログラムは、

コード:

#include "DxLib.h"
#include <math.h>

#define N   8
#define M   3
#define EPS 0.0001

double a[M+1][M+2];
double x[N] = {  -4.0, -2.0, -1.0, 0.0, 1.0, 3.0,  4.0,   6.0 };
double y[N] = { -35.1, 15.1, 15.9, 8.9, 0.1, 0.1, 21.1, 135.0 };

double f(double x)
{
	double y = 0;
	for (int i = 3; i >= 0; i--) y = y * x + a[i][M+1];
	return y;
}

#define X(x) ((int)(O + (x - xa) * W / w))
#define Y(y) ((int)(O + H - (y - ya) * H / h))

void draw()
{
	double dx = 0.2;
	double xa = -5, xb = 7, ya = f(xa), yb = f(xb);
	double w = xb - xa, h = yb - ya;
	double O = 40, W = 560, H = 400;
	int n = (int)(w / dx);

	int x1 = X(0), y1 = Y(0);
	DrawLine(x1, O, x1, O + H, 0x80);
	DrawLine(O, y1, O + W, y1, 0x80);

	x1 = X(xa), y1 = Y(f(xa));
	for (int i = 1; i <= n; i++) {
		double x = xa + i * dx;
		int x2 = X(x), y2 = Y(f(x));
		DrawLine(x1, y1, x2, y2, 0xc00000);
		x1 = x2, y1 = y2;
	}
	for (int i = 0; i < N; i++)
		DrawCircle(X(x[i]), Y(y[i]), 3, 0);

	for (int i = 0; i <= M; i++)
		DrawFormatString(10, 20*i+10, 0, "p%d =%7.3lf\n", i, a[i][M+1]);
}

int Jordan()
{
	for (int i = 0; i <= M; i++) {
		double pivot = a[i][i];
		if (fabs(pivot) < EPS) {
			//printf("ピボットが許容範囲誤差以下\n");
			return 1;
		}
		for (int j = 0; j <= M + 1; j++)
			a[i][j] /= pivot;

		for (int k = 0; k <= M; k++) {
			double delta = a[k][i];
			for (int j = 0; j <= M + 1; j++)
				if (k != i) a[k][j] -= delta * a[i][j];
		}
	}
	return 0;
}

int calc()
{
	for (int i = 0; i <= M; i++)
		for (int j = 0; j <= M + 1; j++) {
			a[j][i] = 0;
			for (int k = 0; k < N; k++)
				a[j][i] += pow(x[k], i + j);
		}

	for (int j = 0; j <= M; j++)
		for (int k = 0; k < N; k++)
			a[j][M + 1] += y[k] * pow(x[k], j);

	return Jordan();
}

int WINAPI WinMain(HINSTANCE, HINSTANCE, LPSTR, int)
{
	SetBackgroundColor(255, 255, 255);
    ChangeWindowMode(TRUE);
	DxLib_Init();
	calc();
    draw();
	WaitKey();
    DxLib_End();
    return 0;
}

回答ありがとうございます。

係数を求めるプログラムと3次関数を求めるプログラムは同じということなのでしょうか？
また、「DxLib を使って図を書くプログラム」の方なのですが、学校では「Visual Studio 用開発者コマンド」を使っており、DxLibのやり方が申し訳ないのですがよくわかりません...。

本当に「DxLib
のやり方」が知りたいの ? それとも「課題を解きたい」の ?

何を使って、グラフを出すか、課題で指定されていませんか ? 後出しは、回答者に嫌われます。

「Visual Studio 用開発者コマンド」としか、書かれていませんか ?
もしくは、前の課題では、「なんとか」を使った、とか、ないですか ?

A.M さんが書きました: ↑
6年前
係数を求めるプログラムと3次関数を求めるプログラムは同じということなのでしょうか？

同じじゃないんですか？

コード:

	for (i = 0; i <= M; i++)
		printf("p%d = %7.3lf\n", i, a[i][M + 1]);

これを次のように書き換えれば同じでしょう。

コード:

	printf("y = %.2f %+.2f x %+.2f x^2 %+.2f x^3",
		a[0][M+1], a[1][M+1], a[2][M+1], a[3][M+1]);

A.M さんが書きました: ↑
6年前
また、「DxLib を使って図を書くプログラム」の方なのですが、学校では「Visual Studio 用開発者コマンド」を使っており、DxLibのやり方が申し訳ないのですがよくわかりません...。

Windows API を使って書き直してみました。
コンパイルできますか？

コード:

#pragma comment(lib, "user32.lib")
#pragma comment(lib, "gdi32.lib")

#define _CRT_SECURE_NO_WARNINGS
#include <windows.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#define N   8
#define M   3
#define EPS 0.0001

double a[M + 1][M + 2];
double x[N] = {  -4.0, -2.0, -1.0, 0.0, 1.0, 3.0,  4.0,   6.0 };
double y[N] = { -35.1, 15.1, 15.9, 8.9, 0.1, 0.1, 21.1, 135.0 };

double f(double x)
{
	double y = 0;
	for (int i = 3; i >= 0; i--) y = y * x + a[i][M+1];
	return y;
}

int Jordan()
{
	for (int i = 0; i <= M; i++) {
		double pivot = a[i][i];
		if (fabs(pivot) < EPS) {
			//printf("ピボットが許容範囲誤差以下\n");
			return 1;
		}
		for (int j = 0; j <= M + 1; j++)
			a[i][j] /= pivot;

		for (int k = 0; k <= M; k++) {
			double delta = a[k][i];
			for (int j = 0; j <= M + 1; j++)
				if (k != i) a[k][j] -= delta * a[i][j];
		}
	}
	return 0;
}

int calc()
{
	for (int i = 0; i <= M; i++)
		for (int j = 0; j <= M + 1; j++) {
			a[j][i] = 0;
			for (int k = 0; k < N; k++)
				a[j][i] += pow(x[k], i + j);
		}

	for (int j = 0; j <= M; j++)
		for (int k = 0; k < N; k++)
			a[j][M + 1] += y[k] * pow(x[k], j);

	return Jordan();
}

#define X(x) ((int)(O + (x - xa) * W / w))
#define Y(y) ((int)(O + H - (y - ya) * H / h))

void onPaint(HWND hwnd)
{
    PAINTSTRUCT ps;
    HDC hdc = BeginPaint(hwnd, &ps);

	RECT rect;
	GetClientRect(hwnd, &rect);
	double O = 40;
	double W = rect.right - rect.left - 2 * O;
	double H = rect.bottom - rect.top - 2 * O;
	double dx = 0.2;
	double xa = -5, xb = 7, ya = f(xa), yb = f(xb);
	double w = xb - xa, h = yb - ya;
	int n = (int)(w / dx);

	int x1 = X(0), y1 = Y(0);
	MoveToEx(hdc, x1, O, NULL), LineTo(hdc, x1, O + H);
	MoveToEx(hdc, O, y1, NULL), LineTo(hdc, O + W, y1);

	HPEN pen = CreatePen(PS_SOLID, 1, RGB(192, 0, 0));
	HPEN pen0 = SelectObject(hdc, pen);
	x1 = X(xa), y1 = Y(f(xa));
	MoveToEx(hdc, x1, y1, NULL);
	for (int i = 1; i <= n; i++) {
		double x = xa + i * dx;
		int x2 = X(x), y2 = Y(f(x));
		LineTo(hdc, x2, y2);
		x1 = x2, y1 = y2;
	}
	SelectObject(hdc, pen0);

	HBRUSH br = GetStockObject(BLACK_BRUSH);
	HBRUSH br0 = SelectObject(hdc, br);
	for (int i = 0; i < N; i++) {
		x1 = X(x[i]), y1 = Y(y[i]);
		Ellipse(hdc, x1-2, y1-2, x1+2, y1+2);
	}
	SelectObject(hdc, br0);

	char buf[256];
	sprintf(buf, "y = %.2f %+.2f x %+.2f x^2 %+.2f x^3",
			a[0][M+1], a[1][M+1], a[2][M+1], a[3][M+1]);
	TextOut(hdc, 10, 10, buf, strlen(buf));

	DeleteObject(pen);
    EndPaint(hwnd, &ps);
}

LRESULT CALLBACK WindowProc(HWND hwnd, UINT msg, WPARAM wp, LPARAM lp)
{
    switch (msg) {
    case WM_PAINT:   onPaint(hwnd); break;
    case WM_DESTROY: PostQuitMessage(0); break;
    default:  return DefWindowProc(hwnd, msg, wp, lp);
    }
    return 0;
}

int WINAPI WinMain(HINSTANCE h, HINSTANCE p, LPSTR cl, int cs)
{
    WNDCLASS wc = { 0, WindowProc, 0, 0, h, NULL,
		LoadCursor(NULL, IDC_ARROW), (HBRUSH)(COLOR_WINDOW+1), NULL, "wc" };
    if (!RegisterClass(&wc)) return FALSE;
    HWND hwnd = CreateWindowEx(WS_EX_COMPOSITED, "wc", cl,
			WS_OVERLAPPEDWINDOW | WS_VISIBLE,
			80, 60, 640, 480, NULL, NULL, h, NULL);
	if (!hwnd) return 0;

	calc();

	ShowWindow(hwnd, cs);
	UpdateWindow(hwnd);

    MSG msg;
    while (GetMessage(&msg, NULL, 0, 0)) {
        TranslateMessage(&msg);
        DispatchMessage(&msg);
    }
    return msg.wParam;
}

ソールファイルの拡張子を .c にしてコンパイルすると
問題ないのですが、.cpp にして C++ としてコンパイルすると、
型変換のエラーになるようです。
SelectObject や GetStockObject の前に
(HPEN) や (HBRUSH) のキャストを入れてください。

コード:

 	HPEN pen0 = (HPEN)SelectObject(hdc, pen);

 	HBRUSH br = (HBRUSH)GetStockObject(BLACK_BRUSH);
 	HBRUSH br0 = (HBRUSH)SelectObject(hdc, br);

> 係数を求めるプログラムと3次関数を求めるプログラムは同じということなのでしょうか？

あなたが求めている「係数」というは，一体何の係数なのでしょうか？

> 最小二乗近似のプログラムを作成し

最小二乗の側は解法とかが何か指定されているんでしょうか？
（あるいは「このライブラリ使って」とか何とか指定がありそうなものだが…
　あと，グラフまで自前プログラムで描画しろってことなんだろうか…？）

ウィンドウのサイズを変更しても再描画されないので、
case WM_SIZE: を追加しました。

コード:

LRESULT CALLBACK WindowProc(HWND hwnd, UINT msg, WPARAM wp, LPARAM lp)
{
    switch (msg) {
    case WM_PAINT: onPaint(hwnd); break;
    case WM_SIZE: InvalidateRect(hwnd, NULL, FALSE); break;
    case WM_DESTROY: PostQuitMessage(0); break;
    default:  return DefWindowProc(hwnd, msg, wp, lp);
    }
    return 0;
}

Windows API を使用せず、コンソールアプリにするなら、

コード:

#include <stdio.h>   // printf
#include <string.h>  // memset
#include <math.h>    // fabs

#define N   8
#define M   3
#define EPS 0.0001

double a[M + 1][M + 2];
double x[N] = {  -4.0, -2.0, -1.0, 0.0, 1.0, 3.0,  4.0,   6.0 };
double y[N] = { -35.1, 15.1, 15.9, 8.9, 0.1, 0.1, 21.1, 135.0 };

double f(double x)
{
	double y = 0;
	for (int i = 3; i >= 0; i--) y = y * x + a[i][M+1];
	return y;
}

int Jordan()
{
	for (int i = 0; i <= M; i++) {
		double pivot = a[i][i];
		if (fabs(pivot) < EPS) {
			printf("ピボットが許容範囲誤差以下\n");
			return 1;
		}
		for (int j = 0; j <= M + 1; j++)
			a[i][j] /= pivot;

		for (int k = 0; k <= M; k++) {
			double delta = a[k][i];
			for (int j = 0; j <= M + 1; j++)
				if (k != i) a[k][j] -= delta * a[i][j];
		}
	}
	return 0;
}

int calc()
{
	for (int i = 0; i <= M; i++)
		for (int j = 0; j <= M + 1; j++) {
			a[j][i] = 0;
			for (int k = 0; k < N; k++)
				a[j][i] += pow(x[k], i + j);
		}

	for (int j = 0; j <= M; j++)
		for (int k = 0; k < N; k++)
			a[j][M + 1] += y[k] * pow(x[k], j);

	return Jordan();
}

#define WIDTH  80
#define HEIGHT 28

char screen[HEIGHT][WIDTH];

void clearScreen(void) { memset(screen, ' ', sizeof screen); }

void drawScreen(void)
{
	for (int i = 0; i < HEIGHT; i++) printf("%.79s\n", screen[i]);
}

void drawChar(unsigned x, unsigned y, char c)
{
	if (x < WIDTH && y < HEIGHT) screen[y][x] = c;
}

void drawString(int x, int y, const char *s)
{
	while (*s) drawChar(x++, y, *s++);
}

void drawLine(int x1, int y1, int x2, int y2, char c)
{
	int t;
	if (x1 == x2) {
		if (y2 < y1) t = y2, y2 = y1, y1 = t;
		for (; y1 <= y2; y1++) drawChar(x1, y1, c);
	}
	else if (y1 == y2) {
		if (x2 < x1) t = x2, x2 = x1, x1 = t;
		for (; x1 <= x2; x1++) drawChar(x1, y1, c);
	}
	else {
		drawChar(x1, y1, c);
		drawChar(x2, y2, c);
	}
}

#define X(x) ((int)(O + (x - xa) * W / w))
#define Y(y) ((int)(O + H - (y - ya) * H / h))

void draw(void)
{
	clearScreen();

	double O = 5;
	double W = 80 - 2 * O;
	double H = 24;
	double dx = 0.2;
	double xa = -5, xb = 7, ya = f(xa), yb = f(xb);
	double w = xb - xa, h = yb - ya;
	int n = (int)(w / dx);

	int x1 = X(0), y1 = Y(0);
	drawLine(x1, 2, x1, H + 2, '|');
	drawLine(0, y1, O + W, y1, '-');
	drawChar(x1, y1, '+');

	x1 = X(xa), y1 = Y(f(xa));
	for (int i = 1; i <= n; i++) {
		double x = xa + i * dx;
		int x2 = X(x), y2 = Y(f(x));
		drawLine(x1, y1, x2, y2, '*');
		x1 = x2, y1 = y2;
	}
	for (int i = 0; i < N; i++) {
		drawChar(X(x[i]), Y(y[i]), 'O');
	}

	char buf[256];
	sprintf(buf, "y = %.2f %+.2f x %+.2f x^2 %+.2f x^3",
			a[0][M+1], a[1][M+1], a[2][M+1], a[3][M+1]);
	drawString(2, 6, buf);

	drawScreen();
}

int main(void)
{
	calc();
	draw();
}

プログラマ専用SNS ミクプラ

数値計算・近似式の問題がわかりません。

数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。

Re: 数値計算・近似式の問題がわかりません。