// To solve the equations described at // http://cap-lore.com/MathPhys/Implicit/1D.html // In routine solve argument n is the number of unknowns, excludng the ends. typedef double D; typedef struct{D x; D y; D z;} mr; void solve(mr * L, int n, D fh, D lh, D * ans){ void qry(D A, D B, int k){ if(k==n) return; else { D r = 1./(1. - L[k].x*B), a = (L[k].x*A + L[k].y)*r, b=L[k].z*r; qry(a, b, k+1); ans[k] = a + b*ans[k+1];}} qry(fh, 0, 0);} #if 1 #include #include int main(){ D fh = 4., lh = 5.; void tx(mr * t, int n){ D ans[n+2]; ans[0]=fh; ans[n+1]=lh; solve(t, n, fh, lh, ans+1); {int j=n+2; while(j--) { printf("%15.10f", ans[j]); if(j&&j<=n) printf(" %15.10f %15.10f %15.10f", t[j-1].x, t[j-1].y, t[j-1].z); printf("\n");}} {int j = n; while(j--) { printf ("%e ", t[j].x*ans[j] + t[j].y + t[j].z*ans[j+2] - ans[j+1]);}} printf("\n");} {mr t[] = {{1,3,2}, {2,7,3}}; tx(t, 2);} {mr t[] = {{0.03, 3, 0.02}, {0.034, 3.6, 0.28}, {0.022, 4, 0.018}}; tx(t, 3);} return 0;} /* 5.0000000000 -12.0000000000 2.0000000000 7.0000000000 3.0000000000 -17.0000000000 1.0000000000 3.0000000000 2.0000000000 4.0000000000 0.000000e+00 0.000000e+00 5.0000000000 4.1974632083 0.0220000000 4.0000000000 0.0180000000 4.8846912884 0.0340000000 3.6000000000 0.2800000000 3.2176938258 0.0300000000 3.0000000000 0.0200000000 4.0000000000 -8.881784e-16 0.000000e+00 0.000000e+00 */ #endif