Investigate stability of the iterative semi-implicit algorithm, demonstrating that it is symplectic for an arbitrary number of iterations.GWH 6/29/2007This gives the details of the proof at the end of Sec. 2.3 of "An iterative semi-implicit scheme with robust damping", N.F. Loureiro and G.W. Hammett, submitted to the J. of Comp. Phys., http://www.arxiv.org/physics/...QyZJKHJlc3RhcnRHJSpwcm90ZWN0ZWRHIiIiLUkra2VybmVsb3B0c0dGJDYjSSh2ZXJzaW9uRzYlRiRJKF9zeXNsaWJHNiJGK0YlQyQtSSdNYXRyaXhHNiI2IzckNyQiIiIiIiM3JCIiJCIiJUYpStart with the iteration matrix B, defined by Eq.39:
QyQ+SSJCRzYiLUknTWF0cml4RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiM3JDckKipJI2R0R0YlIiIjSSJmR0YlIiIiSSJnR0YlRjItSSIqR0YpNiQiIiUsJkkiTEdGJUYyRjJGMiEiIioqRi9GMkYzRjJGOUYyLUY1NiRGMEY4Rjo3JCooRi9GMkYxRjJGPEY6KiZGOUYyRjhGOkYyQyQmSSJCRzYiNiQiIiIiIiNGJw==QyQtSSV3aXRoRzYiNiNJLkxpbmVhckFsZ2VicmFHNiRJKF9zeXNsaWJHRiUlKnByb3RlY3RlZEchIiI=QyQ+NiRJInZHNiJJIlVHRiYtSS1FaWdlbnZlY3RvcnNHRiY2I0kiQkdGJiIiIg==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzZWLUkjbW9HRiQ2MFEiI0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0Y0LyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGNC8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSVmb3JtR1EhRicvJSdsc3BhY2VHUSQwZW1GJy8lJ3JzcGFjZUdGRi8lKG1pbnNpemVHUSIxRicvJShtYXhzaXplR1EpaW5maW5pdHlGJy1GLDYwUSJ+RidGL0YyRjVGN0Y5RjtGPUY/RkFGREZHRklGTC1JI21pR0YkNiVRInZGJy8lJ2l0YWxpY0dRJXRydWVGJy9GMFEnaXRhbGljRidGTy1GUzYlUSNpc0YnRlZGWUZPLUZTNiVRJHRoZUYnRlZGWUZPLUZTNiVRJWxpc3RGJ0ZWRllGTy1GUzYlUSNvZkYnRlZGWUZPLUZTNiVRLGVpZ2VudmFsdWVzRidGVkZZLUYsNjBRIixGJ0YvRjIvRjZGWEY3RjlGO0Y9Rj8vRkJRJmluZml4RidGRC9GSFEzdmVyeXRoaWNrbWF0aHNwYWNlRidGSUZMRk8tRiw2MFEkYW5kRidGL0YyRjVGN0Y5RjtGPUY/RkFGREZHRklGTEZPRmhuRk8tRlM2JVEoY29sdW1uc0YnRlZGWUZPRl5vRk8tRlM2JVEiVUYnRlZGWUZPLUZTNiVRJGFyZUYnRlZGWUZPRmhuRk8tRlM2JVEuY29ycmVzcG9uZGluZ0YnRlZGWUZPLUZTNiVRLWVpZ2VudmVjdG9yc0YnRlZGWS1JJ21zcGFjZUdGJDYmLyUnaGVpZ2h0R1EmMC4wZXhGJy8lJndpZHRoR1EmMC4wZW1GJy8lJmRlcHRoR0ZjcS8lKmxpbmVicmVha0dRKG5ld2xpbmVGJy1GLDYxRi4vJStleGVjdXRhYmxlR0Y0Ri9GMkY1RjdGOUY7Rj1GP0ZBRkRGR0ZJRkwtRiw2MUZRRl5yRi9GMkY1RjdGOUY7Rj1GP0ZBRkRGR0ZJRkwtRlM2JlEmY2hlY2tGJ0ZWRl5yRllGYHItRlM2JkZqbkZWRl5yRllGYHItRlM2JlEsZWlnZW52ZWN0b3JGJ0ZWRl5yRllGYHItRlM2JlEld2l0aEYnRlZGXnJGWUZgci1GUzYmUSJhRidGVkZeckZZRmByLUZTNiZRJXplcm9GJ0ZWRl5yRllGYHItRlM2JlErZWlnZW52YWx1ZUYnRlZGXnJGWUZecUZcckZgci1JKG1mZW5jZWRHRiQ2JS1GIzY1LUZTNiZRKGNhdXRpb25GJ0ZWRl5yRlktRiw2MVEiOkYnRl5yRi9GMkY1RjdGOUY7Rj1GP0Zoby9GRVEvdGhpY2ttYXRoc3BhY2VGJy9GSEZidEZJRkxGYHJGYHItRlM2JlEmbWFwbGVGJ0ZWRl5yRllGYHItRlM2JlEoY2hhbmdlc0YnRlZGXnJGWUZgckZlckZgci1GUzYmUSZvcmRlckYnRlZGXnJGWUZgci1GUzYmRmBvRlZGXnJGWUZgckZlckZgci1GUzYmRmNvRlZGXnJGWUZgci1GUzYmUSpzb21ldGltZXNGJ0ZWRl5yRllGXnJGLw==QyY+SSx6ZWlnZW5pbmRleEc2IiIiIkYmPkkjZTJHRiUtSSdNYXRyaXhHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2IzckNyMmSSJVR0YlNiRGJkYkNyMmRjI2JCIiI0YkRiY=QyQ+SSNlM0c2Ii1JMGRlbGF5RG90UHJvZHVjdEc2JCUqcHJvdGVjdGVkRy9JK21vZHVsZW5hbWVHRiVJLFR5cGVzZXR0aW5nRzYkRilJKF9zeXNsaWJHRiU2JEkiQkdGJUkjZTJHRiUiIiI=QyQtSSlzaW1wbGlmeUc2IjYjSSNlM0dGJSIiIg==QyQ+SSREaWFHNiItSSdNYXRyaXhHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2IzckNyQmSSJ2R0YlNiMiIiIiIiE3JEYyJkYvNiMiIiNGMQ==QyQ+SSRlMTBHNiIpSSREaWFHRiVJInBHRiUiIiI=QyQtSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSRlMTBHRigiIiI=Since this didn't work to simplify the matrix power, do it directly:QyQ+SSVEaWFQRzYiLUknTWF0cml4RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiM3JDckKSZJInZHRiU2IyIiIkkicEdGJSIiITckRjQpJkYwNiMiIiNGM0YyQyQ+SSVEaWFQRzYiKihGJCIiIkkmcmF0aW9HRiVGJyksJComLCZJIkxHRiUiIiUqKEkjZHRHRiUiIiNJImZHRiVGJ0kiZ0dGJUYnRidGJywmRi1GJ0YnRichIiIjRidGLkkicEdGJUY1Ric=QyRJIlVHNiIiIiI=QyQ+SSVVaW52RzYiKiRJIlVHRiUhIiIiIiI=QyQsJi1JMGRlbGF5RG90UHJvZHVjdEc2JCUqcHJvdGVjdGVkRy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nRzYkRidJKF9zeXNsaWJHRio2JC1GJTYkSSJVR0YqSSREaWFHRipJJVVpbnZHRioiIiJJIkJHRiohIiJGNA==QyQtSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSIlR0YoIiIiThus we have verified that B can be written as U D U-1.
Now calculate Bp:QyQ+SSNCcEc2Ii1JMGRlbGF5RG90UHJvZHVjdEc2JCUqcHJvdGVjdGVkRy9JK21vZHVsZW5hbWVHRiVJLFR5cGVzZXR0aW5nRzYkRilJKF9zeXNsaWJHRiU2JC1GJzYkSSJVR0YlSSVEaWFQR0YlSSVVaW52R0YlIiIiQyQ+SSNCcEc2Ii1JKXNpbXBsaWZ5RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiNGJCIiIg==
Next calculate the matrix corresponding to the exact Crank-Nicolson (CN) solution. The exact CN solution is of the form:L yn+1 = R ynwhich we will rearrange as yn+1 = CN ynwhere CN = L-1 R. The L and R matrices are:
QyQ+SSRlMTlHNiItSSdNYXRyaXhHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2IzckNyQiIiIsJComSSNkdEdGJUYuSSJnR0YlRi4jISIiIiIjNyQsJComSSJmR0YlRi5GMUYuRjNGLkYuQyQ+SSRlMjBHNiItSSdNYXRyaXhHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2IzckNyQiIiIsJComSSNkdEdGJUYuSSJnR0YlRi4jRi4iIiM3JCwkKiZJImZHRiVGLkYxRi5GM0YuRi4=QyQ+SSNDTkc2Ii1JMGRlbGF5RG90UHJvZHVjdEc2JCUqcHJvdGVjdGVkRy9JK21vZHVsZW5hbWVHRiVJLFR5cGVzZXR0aW5nRzYkRilJKF9zeXNsaWJHRiU2JCokSSRlMTlHRiUhIiJJJGUyMEdGJSIiIg==QyQ+SSNDTkc2Ii1JKXNpbXBsaWZ5RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiNGJCIiIg==The matrix corresponding to the solution from the first iteration (calculating p=1), as given by Eq. 34, is:
QyQ+SSNBMEc2Ii1JJ01hdHJpeEc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkdGJTYjNyQ3JCwmIiIiRi8qKkkiZkdGJUYvSSJnR0YlRi9JI2R0R0YlIiIjLCZJIkxHRiVGNEY0Ri8hIiJGLyoqRjNGL0YyRi8sKCooRjNGNEYxRi9GMkYvRi9GNiIiJUY7Ri9GLy1JIipHRik2JEY7LCZGNkYvRi9GL0Y3NyQqKEY/RjdGM0YvRjFGL0YuRi8=Try directly proving the result now, without changing to the basis of B:QyQ+SSNBcEc2Ii1JKXNpbXBsaWZ5RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiMsJkkjQ05HRiUiIiItSTBkZWxheURvdFByb2R1Y3RHNiRGKS9JK21vZHVsZW5hbWVHRiVJLFR5cGVzZXR0aW5nR0YoNiRJI0JwR0YlLCZJI0EwR0YlRi5GLSEiIkYuRi4=QyQ+SSR2QXBHNiItSSxFaWdlbnZhbHVlc0dGJTYjSSNBcEdGJSIiIg==QyQ+SSR2QXBHNiItSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkdGJTYjRiQiIiI=QyQ+SSRlMjlHNiItSSlhc3N1bWluZ0dGJTYkNyMtSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkdGJTYjLUklc3Vic0dGLTYlL0kiZkdGJSwkKiZJJm9tZWdhR0YlIiIjSSJnR0YlISIiRjovSSNkdEdGJSomSSRkdG5HRiUiIiItSSRhYnNHRi02I0Y3RjpJJHZBcEdGJTcjJ0Y3SSVyZWFsR0YsRj8=QyQ+SSRlMzBHNiItSSlhc3N1bWluZ0dJKF9zeXNsaWJHRiU2JDcjLUkpc2ltcGxpZnlHNiQlKnByb3RlY3RlZEdGKDYjSSRlMjlHRiU3JTIiIiFJJGR0bkdGJTJGM0kiTEdGJTIqJEY0IiIjLCRGNiIiJSIiIg==For simplicity, select one of the eigenvalues for simplificiation. Repeat this calculation checking the other eigenvalue. (For this problem, I'm pretty sure the eigenvalues come in complex conjugate pairs, so that showing that one of the eigenvalues satisfies |lambda|=1 is sufficient. But it is easy enough to check both eigenvalues.)
Select each one here:QyQ+SSRlMzFHNiImSSRlMzBHRiU2IyIiIyIiIg==Maple needs more help recognizing that the argument of the square root is always negative (if the semi-implicit operator is sufficiently large), and so the sqrt(-arg) can be written as I*sqrt(arg). Select out just the argument of the square root (depending on the order Maple has put the terms, one needs to change the argument to the "op" command below):
QyQ+SSZlMzFiMUc2Ii1JI29wRyUqcHJvdGVjdGVkRzYkNyYiIiMiIihGKyIiIkkkZTMxR0YlRi0=QyQ+SSZlMzFiMkc2Ii1JKXNpbXBsaWZ5RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiNJJmUzMWIxR0YlIiIiQyQ+SSZlMzFiM0c2Ii1JJ2ZhY3Rvckc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkdGJTYjSSZlMzFiMkdGJSIiIg==It is straightforward to prove that the above (the argument of the square root in the eigenvalue) is negative (that the first factor is always positive and the second factor always negative) for arbitrary 0<ratio<1, as long as 4*L>dtn^2. I.e., a sufficient condition for the argument to be negative is omega_hat^2 > omega^2.However, it is actually a little more general than that. I have been able to show that the condition that the first iteration be symplectic 4 + (omega_hat^2 - omega^2) dt^2 > 0is also a sufficient condition for all of the following iterations to be symplectic. This is a little more complicated to prove, and one must consider the possibility of -1<ratio<0 as a special case, but with a few pages of checking the various cases, one can show that this is true.Now knowing that the argument of the square root is always negative, we will make some substitutions to isolate it and then simplify it:QyQ+SSRlMzJHNiItSSdmYWN0b3JHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2I0kkZTMxR0YlIiIiQyQ+SSZlMzJiMUc2Ii1JJXN1YnNHJSpwcm90ZWN0ZWRHNiQvKiYsKiooSSZyYXRpb0dGJSIiIkkiTEdGJUYvSSRkdG5HRiUiIiMiIiUqJkYuRi9GMUYzISIiIiM7Ri9GMEY2Ri8sKiomRi5GL0YwRi9GMyomRi5GL0YxRjJGNSEiJUYvRjBGOkYvLCRJJGFyZ0dGJUY1SSRlMzJHRiVGLw==QyQ+SSRlMzNHNiItSSlhc3N1bWluZ0dJKF9zeXNsaWJHRiU2JDcjLUkpc2ltcGxpZnlHNiQlKnByb3RlY3RlZEdGKDYjSSZlMzJiMUdGJTcjMiIiIUkkYXJnR0YlIiIiCalculate the squared magnitude of the eigenvalue |lambda|^2 :QyQ+SSRlMzRHNiIqJkkkZTMzR0YlIiIiLUkqY29uanVnYXRlRyUqcHJvdGVjdGVkRzYjRidGKEYoQyQ+SSRlMzVHNiItSSlhc3N1bWluZ0dJKF9zeXNsaWJHRiU2JDcjLUkpc2ltcGxpZnlHNiQlKnByb3RlY3RlZEdGKDYjSSRlMzRHRiU3JidJJGR0bkdGJUklcmVhbEdGLSdJIkxHRiVGNCdJJnJhdGlvR0YlRjQyIiIhSSRhcmdHRiUiIiI=QyQ+SSRlMzZHNiItSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkdGJTYjLUknZXhwYW5kR0YpNiNJJGUzNUdGJSIiIg==QyQ+SSRlMzdHNiItSSVzdWJzRyUqcHJvdGVjdGVkRzYkL0kkYXJnR0YlLCQqJiwqKihJJnJhdGlvR0YlIiIiSSJMR0YlRjFJJGR0bkdGJSIiIyIiJSomRjBGMUYzRjUhIiIiIztGMUYyRjhGMSwqKiZGMEYxRjJGMUY1KiZGMEYxRjNGNEY3ISIlRjFGMkY8RjFGN0kkZTM2R0YlRjE=QyQtSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSRlMzdHRigiIiI=Q.E.D.
This completes the proof that the eigenvalues for an arbitrary number of iterations satisfy |lambda_i| = 1, so the algorithm is symplectic.