Witam, dostałem takie zadanie - mam interpolować węzły zadane w pliku wielomianem lagrange'a stopnia trzeciego, a następnie policzyć całkę na zadanym przedziale metodą simpsona. Póki co udało mi się napisać funkcję interpolującą, ale nie jestem pewien co do jej poprawności. Mianowicie wielomian ma silnie oscylacyjny charakter i nie przypomina profilu lotniczego, który przypominać powinien. Prosiłbym w miarę możliwości o sprawdzenie kodu i uwagi, co jest nie tak. Czy w ogóle dobrze ten wielomian wyznaczyłem? Idea jest taka, aby dla każdego argumentu znaleźć cztery najbliższe węzły i wykorzystać je do interpolacji.
#include <stdio.h>
#include <math.h>
#include "winbgi2.h"
//interpoluje podajac do funkcji lagr cztery najblizsze dla danego argumentu wezly
float lagr(float x[], float y[], int n1, float xx);
float simpson(float a, float b, float(*function) (float x), int n);
void main()
{
int n, n1, i, m1, m2, m3, m4, k, t;
float xp, xk, dx, xb, xx, yy,x1,y1;
//float a,b,ct;
float x[100], y[100], min, min2, min3, min4;
float x4[5], y4[5];
graphics(600, 600);
FILE *plik1 = fopen("profil.dat", "r");
fscanf_s(plik1, "%d", &n);
printf_s("ilosc wezlow interpolacji n=%d\n", n);
for (i = 1; i <= n; i++)
{
fscanf_s(plik1, "%f%f", &x[i], &y[i]);
}
xp = 1; xk = 95; dx = 1;
n1 = (xk - xp) / dx + 1;
x1 = xp; y1 = 1.5;
xx = xp;
scale(.0, .0, 100, 10);
for (i = 1; i <= n1; i++)
{
min = x[n];//najwiekszy argument interpolacji
//szukam najbliższego wezla intepolacji
for (k = 1; k <= n; k++){
if (min>fabs(xx - x[k])){
min = fabs(xx - x[k]);
m1 = k;
}
}
min2 = x[n];
//szukam 2. najbliższego wezla intepolacji
for (k = 1; k <= n; k++){
if ((min2>fabs(xx - x[k]) && (fabs(xx - x[k])>min))){
min2 = fabs(xx - x[k]);
m2 = k;
}
}
min3 = x[n];
//szukam 3. najbliższego wezla intepolacji
for (k = 1; k <= n; k++){
if ((min3>fabs(xx - x[k]) && (fabs(xx - x[k])>min2))){
min3 = fabs(xx - x[k]);
m3 = k;
}
}
min4 = x[n];
//szukam 4. najbliższego wezla intepolacji
for (k = 1; k <= n; k++){
if ((min4>fabs(xx - x[k]) && (fabs(xx - x[k])>min3))){
min4 = fabs(xx - x[k]);
m4 = k;
}
}
printf("\nx=%f x1n=%f x2n=%f", xx, x[m1], x[m2]);
x4[1] = x[m1];
x4[4] = x[m4];
y4[1] = y[m1];
y4[4] = y[m4];
x4[2] = x[m2];
y4[2] = y[m3];
x4[3] = x[m3];
y4[3] = y[m3];
yy = lagr(x4, y4, 4, xx);
line(x1, y1, xx, yy);
x1 = xx; y1 = yy;
xx = xx + dx;
}
fclose(plik1);
wait();
}