## Program Code for Solution (Inversion) of Linear System of Equations/****************************************************************************** * This C/C++ function solves a linear system of equations * * ΣkAi,kxk=yi * * for the unknown vector xk, with yi the given vector (the data) and Ai,k * the given coefficient matrix. * This is achieved by algebraically solving the equation for the component xi, i.e. * * xi=(yi -Σk(≠i)Ai,kxk)/Ai,i * * and then iterating this equation using an initial guess for xk. However, in order to * stabilize the algorithm in case of ill-conditioned matrices Ai,k an additional vector ai.xi * is added to the original equation, which thus results in the iteration equation * * xi(n)=(yi +ai.xi(n-1) -Σk(≠i)Ai,kxk(n-1))/(Ai,i+ai) . * * It was found that convergence of the iteration is achieved even for ill-conditioned cases * if ai is chosen such that Ai,i+ai results in a strictly diagonally dominant matrix. * The best way (i.e. providing the best compromise between accuracy of the retrieval and number * of iterations) proved indeed to be to set ai equal to the sum over the absolute values * of the off-diagonal row elements, i.e. * * ai= Σk(≠i)|Ai,k| * * * Input parameters are: * N: number of vector elements, and dimension of the coefficient array * y: vector of data points (with dimension N) * A1: pointer to first element of the coefficient array * epsi: desired relative accuracy (should be chosen not smaller than * the relative noise in the data y) * * Output parameters are: * x: solution vector (with dimension N) * ni: number of iterations * * * The function is called in the form inv(N,y,&A[0][0],epsi,x,ni); * * * * Note that the function not only checks whether the solution vector x has converged, * but also whether the resulting model data have converged to the original data y (the * latter should actually be considered the primary criterion, as there may not be a * unique solution for x considering the given noise level (the solution should be * in particular considered with caution in this sense if convergence was achieved * already after 1 iteration)). * In case the algorithm does not converge, it might help to rearrange the order of * the indexing for the matrix Ai,k (e.g. to AN-i,k). In any case, the matrix should * be reasonably well behaved (random coefficients are unlikely to lead to a * convergence of the algorithm; also note that I have not tested at all whether * the algorithm actually works also in case of mixed positive and negative coefficients * (even though I have formally taken this case into account for the conditioning of the * iteration matrix)). * * This is a rather new algorithm I have developed, and I have tested this so far only * in context with a radiative transfer inversion problem (where it has worked * reasonably well; see http://www.plasmaphysics.org.uk/retrieval.htm ). I would * therefore appreciate any feedback regarding its suitability in other contexts * (for contact details see http://www.plasmaphysics.org.uk/feedback.htm ). * * * Thomas Smid, April 2008 *********************************************************************************/ void inv(int N, double y[], double *A1, double epsi, double x[], int& ni) { int i,k; double A[N][N],a[N],xx[N],yy[N],convgc,dev,devmax; for (i=0; i<=N-1; i=i+1) { for (k=0; k<=N-1; k=k+1) { A[i][k]=A1[N*i+k]; } } for (i=0; i<=N-1; i=i+1) { a[i]=0; for (k=0; k<=N-1; k=k+1) { if (k!=i) {a[i]=a[i]+fabs(A[i][k]);} } if (A[i][i]<=0) {a[i]=1.1*a[i]+fabs(A[i][i]);} A[i][i]=A[i][i]+a[i]; x[i]=y[i]/A[i][i]; } convgc=1; ni=0; while (convgc >= epsi) { for (i=0; i<=N-1; i=i+1) { xx[i]=x[i]; } for (i=0; i<=N-1; i=i+1) { x[i]=0; for (k=0; k<=N-1; k=k+1) { if (k!=i) { x[i]=x[i] +xx[k]*A[i][k]; } } x[i]=(y[i]+a[i]*xx[i]-x[i])/A[i][i] ; } for (i=0; i<=N-1; i=i+1) { yy[i]=0; for (k=0; k<=N-1; k=k+1) { yy[i]=yy[i]+ x[k]*A[i][k]; } yy[i]=yy[i]- x[i]*a[i]; } devmax=0; for (i=1; i<=N-1; i=i+1) { if (xx[i]!=0) { dev=fabs(x[i]/xx[i]-1.0); if (dev>devmax) devmax=dev; } if (y[i]!=0) { dev=fabs(yy[i]/y[i]-1.0); if (dev>devmax) devmax=dev; } } convgc=devmax; ni=ni+1; } } |

Programs Home |

Home |