Image:Lagrangepolys.png
From Wikipedia, the free encyclopedia
Lagrangepolys.png (17KB, MIME type: image/png)
Lagrange interpolation polynomials
Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
Subject to disclaimers.
//
This work is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
This work is distributed in the hope that it will be useful, but without any warranty; without even the implied warranty of merchantability or fitness for a particular purpose. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825
#define PI2 (PI*2)
#define SQ2 1.414213562373095048801688724209698078569671875376948073176679737990732478462
#define FI 1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204
#define SX 201
#define SY 201
#define BPL ((SX*3+3)&~3)
unsigned char bhdr[54]={
0x42, 0x4D, 0x36, 0x00, 0x30, 0x00, 0x00, 0x00, 0x00, 0x00, 0x36, 0x00, 0x00, 0x00, 0x28, 0x00,
0x00, 0x00, 0x00, 0x04, 0x00, 0x00, 0x00, 0x04, 0x00, 0x00, 0x01, 0x00, 0x18, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x30, 0x00, 0x12, 0x0B, 0x00, 0x00, 0x12, 0x0B, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
unsigned char po[BPL];
double gr[SY][SX][3];
void drawit();
int main(int a, char **b) {
FILE *o;
int x, y, c;
double t;
char *p;
srand(time(0));
drawit();
p=bhdr+2; *p++=x=54+BPL*SY; *p++=x>>=8; *p++=x>>=8; *p=x>>=8;
p=bhdr+18; *p++=x=SX; *p++=x>>=8; *p++=x>>=8; *p++=x>>=8;
*p++=x=SY; *p++=x>>=8; *p++=x>>=8; *p=x>>=8;
if(!(o=fopen("func.bmp", "wb"))) {
fclose(o);
printf("Couldn't open output file.\n");
return(0);
}
fwrite(bhdr, 54, 1, o);
for(x=SX*3;x<BPL;++x) po[x]=0;
for(y=SY-1;~y;--y) {
for(x=0,p=po;x<SX;++x) for(c=2;~c;--c) *p++=(t=gr[y][x][c])<=0?0:(t>=1?255:t*255);
fwrite(po, BPL, 1, o);
}
fclose(o);
return(0);
}
double minx, miny, maxx, maxy, dlx, dly;
double cr, cg, cb;
double lagx[4]={-9, -4, -1, 7}, lagy[4]={5, 2, -2, 9}; //Points for lagrange polynomial
void func(int f, double x, double *y, double *dydx) {
double x2, x3, x4, x5, x6, x7, z, w;
int i, j;
x2=x*x;
x3=x2*x;
x4=x3*x;
x5=x4*x;
x6=x5*x;
x7=x6*x;
switch(f) {
//x
case 0: *y=x; *dydx=1; break;
//sinh(x)
case 1: *y=(exp(x)-exp(-x))/2; *dydx=(exp(x)+exp(-x))/2; break;
//cosh(x)
case 2: *y=(exp(x)+exp(-x))/2; *dydx=(exp(x)-exp(-x))/2; break;
//tanh(x)
case 3: *y=(exp(x)-exp(-x))/(exp(x)+exp(-x)); *dydx=1-*y**y; break;
//Chebyshev 0-7
case 100: *y=1; *dydx=0; break;
case 101: *y=x; *dydx=1; break;
case 102: *y=2*x2-1; *dydx=4*x; break;
case 103: *y=4*x3-3*x; *dydx=12*x2-3; break;
case 104: *y=8*x4-8*x2+1; *dydx=32*x3-16*x; break;
case 105: *y=16*x5-20*x3+5*x; *dydx=80*x4-60*x2+5; break;
case 106: *y=32*x6-48*x4+18*x2-1; *dydx=192*x5-192*x3+38*x; break;
case 107: *y=64*x7-112*x5+56*x3-7*x; *dydx=448*x6-560*x4+168*x2-7; break;
//Weirdishev 0-7
case 200: *y=1; *dydx=0; break;
case 201: *y=2*x; *dydx=2; break;
case 202: *y=4*x2-1; *dydx=8*x; break;
case 203: *y=8*x3-4*x; *dydx=24*x2-4; break;
case 204: *y=16*x4-12*x2+1; *dydx=64*x3-24*x; break;
case 205: *y=32*x5-32*x3+6*x; *dydx=160*x4-96*x2+6; break;
//Lagrange polynomial, composite and basis
case 300: for(i=0,*y=0,*dydx=0;i<4;++i) { func(301+i, x, &z, &w); *y+=z; *dydx+=w; } break;
case 301: case 302: case 303:
case 304: for(i=0,z=1;i<4;++i) if(i!=f-301) z*=(x-lagx[i])/(lagx[f-301]-lagx[i]);
for(i=0,w=0;i<4;++i) if(i!=f-301&&x-lagx[i]!=0) w+=z/(x-lagx[i]);
*y=z*lagy[f-301]; *dydx=w*lagy[f-301]; break;
//for(j=0,z=0;j<4;++j) if(j!=f-301) for(i=0,w*y=z;
default: *y=100; *dydx=0; break;
}
}
void subp(int x, int y, double r, double g, double b) {
if(x>=0&&y>=0&&x<SX&&y<SY) {
gr[y][x][0]-=r; gr[y][x][1]-=g; gr[y][x][2]-=b;
}
}
void drawdot(double x, double y) {
int ix, iy;
double dx, dy, ax, ay;
x=(x-minx)/(maxx-minx)*(SX-1)+.5;
y=(y-maxy)/(miny-maxy)*(SY-1)+.5;
ix=floor(x); dx=x-ix; iy=floor(y); dy=y-iy;
ax=1-dx; ay=1-dy;
subp(ix-1, iy-1, cr*ax*ay*.05, cg*ax*ay*.05, cb*ax*ay*.05);
subp(ix , iy-1, cr *ay*.05, cg *ay*.05, cb *ay*.05);
subp(ix+1, iy-1, cr*dx*ay*.05, cg*dx*ay*.05, cb*dx*ay*.05);
subp(ix-1, iy , cr*ax *.05, cg*ax *.05, cb*ax *.05);
subp(ix , iy , cr *.05, cg *.05, cb *.05);
subp(ix+1, iy , cr*dx *.05, cg*dx *.05, cb*dx *.05);
subp(ix-1, iy+1, cr*ax*dy*.05, cg*ax*dy*.05, cb*ax*dy*.05);
subp(ix , iy+1, cr *dy*.05, cg *dy*.05, cb *dy*.05);
subp(ix+1, iy+1, cr*dx*dy*.05, cg*dx*dy*.05, cb*dx*dy*.05);
}
void drawhorz(double y) {
int ix, iy;
double dy, ay;
y=(y-maxy)/(miny-maxy)*(SY-1)+.5;
iy=floor(y); dy=y-iy;
ay=1-dy;
for(ix=0;ix<SX;++ix) {
subp(ix , iy-1, cr *ay , cg *ay , cb *ay );
subp(ix , iy , cr , cg , cb );
subp(ix , iy+1, cr *dy , cg *dy , cb *dy );
}
}
void drawvert(double x) {
int ix, iy;
double dx, ax;
x=(x-minx)/(maxx-minx)*(SX-1)+.5;
ix=floor(x); dx=x-ix;
ax=1-dx;
for(iy=0;iy<SY;++iy) {
subp(ix-1, iy , cr*ax , cg*ax , cb*ax );
subp(ix , iy , cr , cg , cb );
subp(ix+1, iy , cr*dx , cg*dx , cb*dx );
}
}
void drawaxes() {
drawhorz(0);
drawvert(0);
}
void drawgrid() {
int a, b;
for(a=ceil(miny/dly)-1;a<=floor(maxy/dly)+1;++a) drawhorz(a*dly);
for(a=ceil(minx/dlx)-1;a<=floor(maxx/dlx)+1;++a) drawvert(a*dlx);
}
void drawfunc(int f) {
double x, y, dydx, pfx, pfy;
pfx=(maxx-minx)/(SX-1);
pfy=(maxy-miny)/(SY-1);
for(x=minx;x<maxx;x+=.1*pfx/sqrt(1+dydx*dydx/((pfy*pfy)/(pfx*pfx)))) {
func(f, x, &y, &dydx);
drawdot(x, y);
}
}
void drawcirc(double x, double y, double r) {
double a, pfx, pfy;
pfx=(maxx-minx)/(SX-1);
pfy=(maxy-miny)/(SY-1);
for(a=0;a<PI2;a+=.1/sqrt(sin(a)*sin(a)/(pfx*pfx)+cos(a)*cos(a)/(pfy*pfy))/r)
drawdot(x+r*cos(a), y+r*sin(a));
}
void drawit() {
int x, y, c;
for(y=0;y<SY;++y) for(x=0;x<SY;++x) for(c=0;c<3;++c) gr[y][x][c]=1;
//Lagrange polynomials, composite and basis (*y)
minx=miny=-10.; maxx=maxy=10.; dlx=dly=1;
cr=.6; cg=.6; cb=.6;
drawfunc(300);
cr=.1; cg=.8; cb=.8;
drawfunc(301); drawcirc(lagx[0], lagy[0], .6);
cr=.8; cg=.8; cb=.1;
drawfunc(302); drawcirc(lagx[1], lagy[1], .6);
cr=.8; cg=.1; cb=.8;
drawfunc(303); drawcirc(lagx[2], lagy[2], .6);
cr=.1; cg=.1; cb=.8;
drawfunc(304); drawcirc(lagx[3], lagy[3], .6);
cr=cg=cb=.8; drawaxes();
cr=cg=cb=.1; drawgrid();
/* //Chebyshev 0, 1, 2, 3, 4, 5//, 6, 7
minx=miny=-5/4.; maxx=maxy=5/4.; dlx=dly=1;
cr=.6; cg=.6; cb=.6;
drawfunc(100);
cr=.1; cg=.8; cb=.8;
drawfunc(101);
cr=.8; cg=.8; cb=.1;
drawfunc(102);
cr=.8; cg=.1; cb=.8;
drawfunc(103);
cr=.1; cg=.1; cb=.8;
drawfunc(104);
cr=.1; cg=.1; cb=.1;
drawfunc(105);
cr=.1; cg=.8; cb=.1;
//drawfunc(106);
cr=.8; cg=.1; cb=.1;
//drawfunc(107);
cr=cg=cb=.8; drawaxes();
cr=cg=cb=.1; drawgrid();
cr=cg=cb=.025; dlx=dly=.1; drawgrid();
*/
/* //Chebyshev polynomials of the weird kind 0, 1, 2, 3, 4, 5
minx=miny=-5/4.; maxx=maxy=5/4.; dlx=dly=1;
cr=.6; cg=.6; cb=.6;
drawfunc(200);
cr=.1; cg=.8; cb=.8;
drawfunc(201);
cr=.8; cg=.8; cb=.1;
drawfunc(202);
cr=.8; cg=.1; cb=.8;
drawfunc(203);
cr=.1; cg=.1; cb=.8;
drawfunc(204);
cr=.1; cg=.1; cb=.1;
drawfunc(205);
cr=cg=cb=.8; drawaxes();
cr=cg=cb=.1; drawgrid();
cr=cg=cb=.025; dlx=dly=.1; drawgrid();
*/
/* //sinh, cosh, tanh
minx=miny=-5; maxx=maxy=5; dlx=dly=1;
cr=.1; cg=cb=.8;
drawfunc(1);
cg=.1; cb=cr=.8;
drawfunc(2);
cb=.1; cr=cg=.8;
drawfunc(3);
cr=cg=cb=.8; drawaxes();
cr=cg=cb=.1; drawgrid();
*/
}
//
File history
Legend: (cur) = this is the current file, (del) = delete
this old version, (rev) = revert to this old version.
Click on date to download the file or see the image uploaded on that date.
- (del) (cur) 23:13, 20 September 2004 . . Cyp (Talk | contribs) . . 201×201 (16,947 bytes) (Lagrange interpolation polynomials {{GFDL}})
- Edit this file using an external application
See the setup instructions for more information.
File links
The following pages on the English Wikipedia link to this file (pages on other projects are not listed):
