c17-3 (779603), страница 3
Текст из файла (страница 3)
The parameters conv, slowc, scalv relate to convergence. Eachinversion of the matrix produces corrections for ne variables at m mesh points. We want theseto become vanishingly small as the iterations proceed, but we must define a measure for thesize of corrections. This error “norm” is very problem specific, so the user might wish torewrite this section of the code as appropriate. In the program below we compute a value forthe average correction err by summing the absolute value of all corrections, weighted by ascale factor appropriate to each type of variable:768Chapter 17.Two Point Boundary Value Problemsfraction of the corrections found by matrix inversion:Y [j][k] → Y [j][k] +slowc∆Y [j][k]max(slowc,err)(17.3.14)void difeq(int k, int k1, int k2, int jsf, int is1, int isf,int indexv[], int ne, float **s, float **y);The only information passed from difeq to solvde is the matrix of derivativess[1..ne][1..2*ne+1]; all other arguments are input to difeq and should not be altered.k indicates the current mesh point, or block number.
k1,k2 label the first and last point inthe mesh. If k=k1 or k>k2, the block involves the boundary conditions at the first or finalpoints; otherwise the block acts on FDEs coupling variables at points k-1, k.The convention on storing information into the array s[i][j] follows that used inequations (17.3.8), (17.3.10), and (17.3.12): Rows i label equations, columns j refer toderivatives with respect to dependent variables in the solution. Recall that each equation willdepend on the ne dependent variables at either one or two points.
Thus, j runs from 1 toeither ne or 2*ne. The column ordering for dependent variables at each point must agreewith the list supplied in indexv[j]. Thus, for a block not at a boundary, the first columnmultiplies ∆Y (l=indexv[1],k-1), and the column ne+1 multiplies ∆Y (l=indexv[1],k).is1,isf give the numbers of the starting and final rows that need to be filled in the s matrixfor this block. jsf labels the column in which the difference equations Ej,k of equations(17.3.3)–(17.3.5) are stored. Thus, −s[i][jsf] is the vector on the right-hand side of thematrix. The reason for the minus sign is that difeq supplies the actual difference equation,Ej,k , not its negative.
Note that solvde supplies a value for jsf such that the differenceequation is put in the column just after all derivatives in the s matrix. Thus, difeq expects tofind values entered into s[i][j] for rows is1 ≤ i ≤ isf and 1 ≤ j ≤ jsf.Finally, s[1..nsi][1..nsj] and y[1..nyj][1..nyk] supply difeq with storagefor s and the solution variables y for this iteration. An example of how to use this routineis given in the next section.Many ideas in the following code are due to Eggleton [1].#include <stdio.h>#include <math.h>#include "nrutil.h"void solvde(int itmax, float conv, float slowc, float scalv[], int indexv[],int ne, int nb, int m, float **y, float ***c, float **s)Driver routine for solution of two point boundary value problems by relaxation. itmax is themaximum number of iterations. conv is the convergence criterion (see text).
slowc controlsthe fraction of corrections actually used after each iteration. scalv[1..ne] contains typicalsizes for each dependent variable, used to weight errors. indexv[1..ne] lists the columnordering of variables used to construct the matrix s[1..ne][1..2*ne+1] of derivatives. (Thenb boundary conditions at the first mesh point must contain some dependence on the first nbvariables listed in indexv.) The problem involves ne equations for ne adjustable dependentvariables at each point.
At the first mesh point there are nb boundary conditions. There are atotal of m mesh points. y[1..ne][1..m] is the two-dimensional array that contains the initialguess for all the dependent variables at each mesh point. On each iteration, it is updated by theSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Thus, when err>slowc only a fraction of the corrections are used, but when err≤slowcthe entire correction gets applied.The call statement also supplies solvde with the array y[1..nyj][1..nyk] containing the initial trial solution, and workspace arrays c[1..ne][1..ne-nb+1][1..m+1],s[1..ne][1..2*ne+1]. The array c is the blockbuster: It stores the unreduced elementsof the matrix built up for the backsubstitution step. If there are m mesh points, then therewill be m+1 blocks, each requiring ne rows and ne-nb+1 columns. Although large, thisis small compared with (ne×m)2 elements required for the whole matrix if we did notbreak it into blocks.We now describe the workings of the user-supplied function difeq.
The synopsis ofthe function is17.3 Relaxation Methods769kmax=ivector(1,ne);ermax=vector(1,ne);k1=1;Set up row and column markers.k2=m;nvars=ne*m;j1=1;j2=nb;j3=nb+1;j4=ne;j5=j4+j1;j6=j4+j2;j7=j4+j3;j8=j4+j4;j9=j8+j1;ic1=1;ic2=ne-nb;ic3=ic2+1;ic4=ne;jc1=1;jcf=ic3;for (it=1;it<=itmax;it++) {Primary iteration loop.k=k1;Boundary conditions at first point.difeq(k,k1,k2,j9,ic3,ic4,indexv,ne,s,y);pinvs(ic3,ic4,j5,j9,jc1,k1,c,s);for (k=k1+1;k<=k2;k++) {Finite difference equations at all point pairs.kp=k-1;difeq(k,k1,k2,j9,ic1,ic4,indexv,ne,s,y);red(ic1,ic4,j1,j2,j3,j4,j9,ic3,jc1,jcf,kp,c,s);pinvs(ic1,ic4,j3,j9,jc1,k,c,s);}k=k2+1;Final boundary conditions.difeq(k,k1,k2,j9,ic1,ic2,indexv,ne,s,y);red(ic1,ic2,j5,j6,j7,j8,j9,ic3,jc1,jcf,k2,c,s);pinvs(ic1,ic2,j7,j9,jcf,k2+1,c,s);bksub(ne,nb,jcf,k1,k2,c);Backsubstitution.err=0.0;for (j=1;j<=ne;j++) {Convergence check, accumulate average erjv=indexv[j];ror.errj=vmax=0.0;km=0;for (k=k1;k<=k2;k++) {Find point with largest error, for each devz=fabs(c[jv][1][k]);pendent variable.if (vz > vmax) {vmax=vz;km=k;}errj += vz;}err += errj/scalv[j];Note weighting for each dependent variable.ermax[j]=c[jv][1][km]/scalv[j];Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).calculated correction. The arrays c[1..ne][1..ne-nb+1][1..m+1] and s supply dummystorage used by the relaxation code.{void bksub(int ne, int nb, int jf, int k1, int k2, float ***c);void difeq(int k, int k1, int k2, int jsf, int is1, int isf,int indexv[], int ne, float **s, float **y);void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k,float ***c, float **s);void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf,int ic1, int jc1, int jcf, int kc, float ***c, float **s);int ic1,ic2,ic3,ic4,it,j,j1,j2,j3,j4,j5,j6,j7,j8,j9;int jc1,jcf,jv,k,k1,k2,km,kp,nvars,*kmax;float err,errj,fac,vmax,vz,*ermax;770Chapter 17.Two Point Boundary Value Problems}nrerror("Too many iterations in solvde");Convergence failed.}void bksub(int ne, int nb, int jf, int k1, int k2, float ***c)Backsubstitution, used internally by solvde.{int nbf,im,kp,k,j,i;float xx;nbf=ne-nb;im=1;for (k=k2;k>=k1;k--) {Use recurrence relations to eliminate remaining deif (k == k1) im=nbf+1;pendences.kp=k+1;for (j=1;j<=nbf;j++) {xx=c[j][jf][kp];for (i=im;i<=ne;i++)c[i][jf][k] -= c[i][j][k]*xx;}}for (k=k1;k<=k2;k++) {Reorder corrections to be in column 1.kp=k+1;for (i=1;i<=nb;i++) c[i][1][k]=c[i+nbf][jf][k];for (i=1;i<=nbf;i++) c[i+nb][1][k]=c[i][jf][kp];}}#include <math.h>#include "nrutil.h"void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k, float ***c,float **s)Diagonalize the square subsection of the s matrix, and store the recursion coefficients in c;used internally by solvde.{int js1,jpiv,jp,je2,jcoff,j,irow,ipiv,id,icoff,i,*indxr;float pivinv,piv,dum,big,*pscl;indxr=ivector(ie1,ie2);pscl=vector(ie1,ie2);je2=je1+ie2-ie1;js1=je2+1;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).kmax[j]=km;}err /= nvars;fac=(err > slowc ? slowc/err : 1.0);Reduce correction applied when error is large.for (j=1;j<=ne;j++) {Apply corrections.jv=indexv[j];for (k=k1;k<=k2;k++)y[j][k] -= fac*c[jv][1][k];}printf("\n%8s %9s %9s\n","Iter.","Error","FAC");Summary of correctionsprintf("%6d %12.6f %11.6f\n",it,err,fac);for this step.if (err < conv) {Point with largest error for each variable canfree_vector(ermax,1,ne);be monitored by writing out kmax andfree_ivector(kmax,1,ne);ermax.return;}17.3 Relaxation Methods771}void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf,int ic1, int jc1, int jcf, int kc, float ***c, float **s)Reduce columns jz1-jz2 of the s matrix, using previous results as stored in the c matrix.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
