01: #include <stdio.h> 02: #include <math.h> 03: #define MAXN 1010 04: 05: void mk_graph (int n, double x[], double x1, double x2, 06: double y[], double y1, double y2); 07: 08: /*=============================================================*/ 09: /* main function */ 10: /*=============================================================*/ 11: int main(void){ 12: double xmin, xmax; 13: int i, nplot; 14: double px[MAXN], py[MAXN]; 15: double x0, x1, x2, y0, y1, y2; 16: 17: xmin=-1.0; 18: xmax=5.0; 19: 20: x0=1; y0=1; 21: x1=3; y1=2; 22: x2=4; y2=5; 23: 24: nplot=1000; 25: for(i=0; i<=nplot; i++){ 26: px[i]=xmin+i*(xmax-xmin)/(nplot); 27: py[i]=(px[i]-x1)*(px[i]-x2)/((x0-x1)*(x0-x2))*y0+ 28: (px[i]-x0)*(px[i]-x2)/((x1-x0)*(x1-x2))*y1+ 29: (px[i]-x0)*(px[i]-x1)/((x2-x0)*(x2-x1))*y2; 30: } 31: 32: mk_graph (nplot, px, xmin, xmax, py, 0.0, 6.0); 33: 34: return 0; 35: 36: } 37: 38: /*==========================================================*/ 39: /* make a graph */ 40: /*==========================================================*/ 41: void mk_graph(int n, double x[], double x1, double x2, 42: double y[], double y1, double y2){ 43: int i; 44: char *data_file; 45: FILE *out; 46: FILE *gp; 47: 48: /* == make a data file ========= */ 49: 50: data_file="gpdata.txt"; 51: 52: out = fopen(data_file, "w"); 53: 54: for(i=0; i<=n; i++){ 55: fprintf(out, "%e\t%e\n", x[i], y[i]); 56: } 57: 58: fclose(out); 59: 60: /* == flowing lines make a graph by using gnuplot == */ 61: 62: gp = popen("gnuplot -persist","w"); 63: 64: fprintf(gp, "reset\n"); 65: fprintf(gp, "set terminal postscript eps color\n"); 66: fprintf(gp, "set output \"graph.eps\"\n"); 67: fprintf(gp, "set grid\n"); 68: 69: /* ------- set x axis ---------*/ 70: 71: fprintf(gp, "set xlabel \"%s\"\n", "x"); 72: fprintf(gp, "set nologscale x\n"); 73: fprintf(gp, "set xrange[%e:%e]\n", x1, x2); 74: 75: /* ------- set y axis ---------*/ 76: 77: fprintf(gp, "set ylabel \"%s\"\n", "y"); 78: fprintf(gp, "set nologscale y\n"); 79: fprintf(gp, "set yrange[%e:%e]\n", y1, y2); 80: 81: /* ------- plat graphs ---------*/ 82: 83: fprintf(gp, "plot \"%s\" using 1:2 with line\n", data_file); 84: 85: fprintf(gp, "set terminal x11\n"); 86: fprintf(gp, "replot\n"); 87: 88: pclose(gp); 89: }
001: #include <stdio.h> 002: #include <math.h> 003: #define MAXN 1010 004: 005: void spline_cal_u(int n, double x[], double y[], double u[]); 006: double spline(int n, double x[], double y[], double u[], double xx); 007: int gauss_jordan(int n, double a[][MAXN], double b[]); 008: void mk_graph(int n, double x[], double x1, double x2, 009: double y[], double y1, double y2); 010: 011: /*=============================================================*/ 012: /* main function */ 013: /*=============================================================*/ 014: 015: int main(){ 016: double x[MAXN], y[MAXN], u[MAXN]; 017: double xmin, xmax; 018: int n,i, nplot; 019: double px[MAXN], py[MAXN]; 020: 021: xmin = -1.0; 022: xmax = 1.0; 023: n=11; 024: 025: for(i=0; i<n; i++){ 026: x[i]=xmin+i*(xmax-xmin)/(n-1); 027: y[i]=1.0/(1.0+25*x[i]*x[i]); 028: } 029: 030: spline_cal_u(n,x,y,u); 031: 032: nplot = 1000; 033: for(i=0; i<=nplot; i++){ 034: px[i]=xmin+i*(xmax-xmin)/(nplot); 035: py[i]=spline(n, x, y, u, px[i]); 036: } 037: 038: mk_graph(nplot, px, xmin, xmax, py, -0.5, 1.5); 039: 040: return 0; 041: 042: } 043: /*=============================================================*/ 044: /* calculation of u[i] */ 045: /*=============================================================*/ 046: void spline_cal_u(int n, double x[], double y[], double u[]){ 047: double a[MAXN][MAXN]; 048: int i, j; 049: 050: for(i=0; i<=n-1; i++){ 051: for(j=0; j<=n-1; j++){ 052: a[i][j]=0.0; 053: } 054: } 055: 056: a[1][1]=2.0*(x[2]-x[0]); 057: a[1][2]=x[2]-x[1]; 058: u[1]=6.0*((y[2]-y[1])/(x[2]-x[1])-(y[1]-y[0])/(x[1]-x[0])); 059: 060: for(i=2; i<=n-2; i++){ 061: a[i][i-1]=x[i]-x[i-1]; 062: a[i][i]=2.0*(x[i+1]-x[i-1]); 063: a[i][i+1]=x[i+1]-x[i]; 064: u[i]=6.0*((y[i+1]-y[i])/(x[i+1]-x[i])- 065: (y[i]-y[i-1])/(x[i]-x[i-1])); 066: } 067: 068: a[n-1][n-2]=x[n-1]-x[n-2]; 069: a[n-1][n-1]=2.0*(x[n]-x[n-2]); 070: u[n-1]=6.0*((y[n]-y[n-1])/(x[n]-x[n-1])- 071: (y[n-1]-y[n-2])/(x[n-1]-x[n-2])); 072: 073: gauss_jordan(n-1, a, u); 074: 075: u[0]=0.0; 076: u[n]=0.0; 077: 078: 079: } 080: 081: /*=============================================================*/ 082: /* spline interpolation */ 083: /*=============================================================*/ 084: double spline(int n, double x[], double y[], double u[], 085: double xx){ 086: int lo, hi, k; 087: double a, b, c, d; 088: double h, yy; 089: 090: lo=0; 091: hi=n; 092: 093: while(hi-lo>1){ 094: k=(hi+lo)/2; 095: if(xx < x[k]){ 096: hi=k; 097: }else{ 098: lo=k; 099: } 100: } 101: 102: a=(u[hi]-u[lo])/(6.0*(x[hi]-x[lo])); 103: b=u[lo]/2.0; 104: c=(y[hi]-y[lo])/(x[hi]-x[lo]) 105: -1.0/6.0*(x[hi]-x[lo])*(2*u[lo]+u[hi]); 106: d=y[lo]; 107: 108: h=xx-x[lo]; 109: yy=a*h*h*h+b*h*h+c*h+d; 110: 111: return(yy); 112: 113: } 114: 115: /*=============================================================*/ 116: /* simple Gauss-Jordan method */ 117: /*=============================================================*/ 118: int gauss_jordan(int n, double a[][MAXN], double b[]){ 119: 120: int ipv, i, j; 121: double inv_pivot, temp; 122: 123: for(ipv=1 ; ipv <= n ; ipv++){ 124: 125: inv_pivot = 1.0/a[ipv][ipv]; 126: for(j=1 ; j <= n ; j++){ 127: a[ipv][j] *= inv_pivot; 128: } 129: b[ipv] *= inv_pivot; 130: 131: for(i=1 ; i<=n ; i++){ 132: if(i != ipv){ 133: temp = a[i][ipv]; 134: for(j=1 ; j<=n ; j++){ 135: a[i][j] -= temp*a[ipv][j]; 136: } 137: b[i] -= temp*b[ipv]; 138: } 139: } 140: } 141: return 0; 142: } 143: 144: /*==========================================================*/ 145: /* make a graph */ 146: /*==========================================================*/ 147: void mk_graph(int n, double x[], double x1, double x2, 148: double y[], double y1, double y2){ 149: int i; 150: char *data_file; 151: FILE *out; 152: FILE *gp; 153: 154: /* == make a data file ========= */ 155: 156: data_file="gpdata.txt"; 157: 158: out = fopen(data_file, "w"); 159: 160: for(i=0; i<=n; i++){ 161: fprintf(out, "%e\t%e\n", x[i], y[i]); 162: } 163: 164: fclose(out); 165: 166: /* == flowing lines make a graph by using gnuplot == */ 167: 168: gp = popen("gnuplot -persist","w"); 169: 170: fprintf(gp, "reset\n"); 171: fprintf(gp, "set terminal postscript eps color\n"); 172: fprintf(gp, "set output \"graph.eps\"\n"); 173: fprintf(gp, "set grid\n"); 174: 175: /* ------- set x axis ---------*/ 176: 177: fprintf(gp, "set xlabel \"%s\"\n", "x"); 178: fprintf(gp, "set nologscale x\n"); 179: fprintf(gp, "set xrange[%e:%e]\n", x1, x2); 180: 181: /* ------- set y axis ---------*/ 182: 183: fprintf(gp, "set ylabel \"%s\"\n", "y"); 184: fprintf(gp, "set nologscale y\n"); 185: fprintf(gp, "set yrange[%e:%e]\n", y1, y2); 186: 187: /* ------- plat graphs ---------*/ 188: 189: fprintf(gp, "plot \"%s\" using 1:2 with line\n", data_file); 190: 191: fprintf(gp, "set terminal x11\n"); 192: fprintf(gp, "replot\n"); 193: 194: pclose(gp); 195: } 196:
last update: