Friday, 13 February 2015

Happy Birthday SDO!

It has been a week of highlights!

News was announced of The University of Sheffield's involvement with the Daniel K Inouye Solar Telescope (DKIST). With a four metre diameter primary mirror  "the telescope will be able to pick up unprecedented detail on the surface of the Sun – the equivalent of being able to examine a £1 coin from 100kms away".  Further information about DKIST is on its own site and wikipedia.
The first total eclipse of the Sun to be visible in Northern Europe since that of 11 August 1999 will take place on 20 March 2015. Further new is provided by The British Astronomical Association.

NASA have prepared a news item about SDO's 5th year - "New Videos Highlight NASA SDO's Fifth Anniversary".

A five year time-lapse video of the sun

On February 11th    launched  NOAA’s New Deep Space Solar Monitoring Satellite on a SpaceX Falcon 9 rocket from Cape Canavarel. An issue with the tracking meant that the launch was delayed from February 8th. The Deep Space Climate Observatory, or DSCOVR, is a spacecraft which will orbit between Earth and the sun, observing and providing advanced warning of particles and magnetic fields emitted by the sun (known as the solar wind) which can affect power grids, communications systems, and satellites close to Earth. From its post at the Lagrange point 1 (or L1), approximately one million miles from Earth.

Tuesday, 10 February 2015

Studying Imagery from the Solar Dynamics Observatory with The Matlab Image Processing Toolbox


In this posting we present an overview if image processing techniques  using the matlab image processing toolbox. We apply these skills to a study of imagery generated by the NASA Solar Dynamics Observatory (SDO). SDO has provided for 5 years now high resolution imagery of our nearest star with a cadence of around 1 second. The video below illustrates the wonders captured by SDO and some of the wildly different eruptive events which can be produced by the sun. These events include solar flares, events where there is an additional ejection of solar material called a coronal mass ejection (CME). The event below also shows more complex moving structures in association with changes in magnetic field lines that loop up into the Sun's atmosphere, the corona. On July 19, 2012, an eruption occurred on the sun that produced all three. A moderately powerful solar flare exploded on the Sun's lower right hand limb, sending out light and radiation. Next came a CME, which shot off to the right out into space. And then, the Sun treated viewers to one of its dazzling magnetic displays -- a phenomenon known as coronal rain.

 The event studied was an M-9.3 flare (just short of an X class, the largest) from an active region right at the Sun's edge (Mar. 12-13, 2014).  The analysis focuses on a group of loops we  examine the properties of the loops over a sequence of times before and after the flare event.

Image processing Tools

 Solar physicists are provided with an armoury of analytical tools, these include SolarSoft (IDL)SunPY (Python)MyExperiment (see references 15 and 16) and the Heliophysics integrated laboratory  (see references 15 and 16).  An excellent overview of solar coronal loops is given in reference 22.

The Matlab image processing toolbox provides a range of functions for
Much of this functionality can take advantage of the performance of graphical processing units for processing imagery. The documentation provided with Matlab and online is excellent, there are many helpful examples. There are also a lot of tutorials which are referenced below (see ref.  1-4). Once data had been imported (from SDO) the  imread function was used to load data into Matlab, it takes as input a filename and returns a handle to the loaded image, a processed image can then be written using the imwrite function. Example syntax is below

imwrite(trans_image, looprotimfile);

There is a range of tools for displaying and exploring images. The simplest is the imshow function. The image tool is particularly useful for anlysing gray scale images it is started using the imtool command. which takes as input the handle to the image. Example syntax is below


imshowpair is used to visualise the difference between two images. There are many steps which can be taken to prepare our selected imagery for analysis. The simplest operations include the ability to select regions of interest, adjust the position and orientation of the image. It is also necessary to modify the colour and contrast within the image so as to enable a more effective study of the features of interest. The geometric transformation tools we used were the imtranslate, imrotate and imcrop functions. These were used to translate rotate and select a region of interest of the image respectively. Each of the functions takes a handle to an input image and returns a handle to the output image. They also take as input the parameters required to translate, rotae and crop the image. Examples are shown below.

imtranslate(li2, [-translation, 0],'FillValues',255);
imcrop(li2, [ xmin ymin width height]);

Matlab provides a range of tools for image enhancement from the contrast adjustment subset we used the histeq function to enhance the image contrast. Before doing histogram equalisation we converted images to grey scale using the rgb2gray function from the graphics toolbox. Basic examples are shown below;


The imsharpen function is a useful tool for bringing out features.  This uses the unsharp masking technique to enhance the image. This may sharpen any edges in the image a Gaussian lowpass filter can be tuned by specifying its diameter in terms of the number of pixels. 

To examine and make measurements from the the imagery we used the imshow and imtool functions. Although we can measure pixel positions and properties with image show the imtool provides tools for measuring the distance between pixels in the image viewer.

The final tools we used for processing were the plot, line and improfile tools. The plot and line graphics functions were used to mark and highlight features of interest. We used the improfile tool to provide pixel-value cross-sections along given line segments. Examples of using these are given below



Matlab also provides the impixel and imhist functions to provide pixel data and present information on the image characteristics. Although we've briefly described a limited subset further details are in the online documentation and in the tutorials  in reference 1-4.

Coronal Loops

Our test case is that of processing a sequence of images of coronal loops in the solar atmosphere. Coronal loops are structures observed in the lower corona and transition region of the solar atmosphere. One of the earliest recorded observations of such structures was by  Angelo Secchi an Italian priest and astronomer. From observation's of the solar eclipses in Spain in 1860  he  presented sketches of coronal loop structures in his 1870 publication Le Soleil: Exposé des Principales Découvertes Modernes (The Sun: Presentation of the Major Modern Discoveries). He organised a further eclipse viewing trip in 1870. One of his sketches is shown below. Secchi's sketch also provides evidence for spicule structures observed on the solar limb. Further images of the solar corona from eclipse viewings are available from reference 12.

The coronal loops are actually spectacular evidence for the magnetic structures in the corona. Loops can extend upto 700Mm in length these structures contain plasma which is heated and rise along the magnetic field lines before falling back again with speeds of up to 100km/s. Loop temperatures can be between 2-3MK and the densities in the range 10^14-10^15 m^-3. The ends of the loops are located in regions of strong magnetic fields near the edge of active regions. Individual loops my last for about a day but a loop system may exist for a few solar rotations. More detailed studies reveal that these structures range in complexity consisting of many strands which can be braided and twisted. With the advent of more sophisticated solar telescopes, coronagraphs and space borne telescopes images of these structures are readily available for example the opening video highlighted coronal rain captured by SDO.

Space based observations of solar phenomena (using missions such as SDO, TRACE etc..) show frequent catastrophic cooling and evacuation of quiescent solar coronal loops over active regions. Analysis of this imagery can show plasma from a temperature of a few million degrees down to less than 100 000 K. Radiative emission from the cooling loops develops initially near the loop tops; later, that cool plasma usually slides down on both sides of the loop.  The evacuation of a coronal loop often occurs after plasma high in the corona has cooled to transition-region or even chromospheric temperatures. The relatively cool material often forms clumps that move at speeds of up to 100 km s−1. The downward acceleration is no more than 80 m s−2, less than a third of the surface gravity (ref. 21). Coronal loops are predominantly heated at the footpoints but their  catastrophic cooling can lead to localized brightenings. Loop heating can lead to an evaporation of plasma into the coronal loop which then cools rapidly due to a loss of thermal equilibrium. The confined region of "condensed'' plasma subsequently falls down under the effect of gravity in the form of a cool plasma blob (ref .17). Thermal nonequilibrium may nonetheless play an important role in prominences and catastrophic cooling events (e.g., coronal rain) that occupy a small fraction of the coronal volume ( ref. 20).


Our study focuses on an M2.5 class flare event which occurred on 12th March 2014 (see NOAA), the event was recorded in the SDO event gallery. The M2.5   flare occurred in NOAA group number 11996  located at N14W76. The region we study is NOAA group 12003 located at N06W37, no flares were reported for this group on this day. The Atmospheric Imaging Assembly (AIA) aboard the Solar Dynamics Observatory has state-of-the-art spatial resolution and shows the most detailed images of coronal loops ever observed. The series of coronal filters peak at different temperatures, which span the range of active regions. These features represent a significant improvement over earlier coronal imagers. There are many sources for gathering and exploring solar data including.

SDO AIA/HMI data access
Virtual solar observatory (VSO)

For our study we used the SDO data browser. We gathered two sets of data corresponding to the AIA193 telescope and the AIA171 telescope and including the output from the Helioseismic and magnetic imager. For the AIA193 data set we collected 184 images starting at 00:40:07 on 12th March 2014 and finishing at 22:40:07 in 13th March 2014. For the AIA171 (with HMI data) we collected 121 images starting at 00:29:42 on 12th March 2014 and finishing at 22:43:11 on 13th march 2014. For both sets we collected images with a resolution of 4096x4096 pixels. It would have been helpful to gather data using the VSO and to gather FITS data froom SDO. Using the flexible image transport system (FITS) tool box available with Matlab it would be straightforward to extract and  use SDO metadata in our analysis.  The event may also be viewed using helioviewer (2014-03-12T22:36:37) this was an easy way to view the metadata from SDO.

A video of the region we studied generated using helioviewer is shown below for the AIA193 detector.

The videos below were generated from the sequence of images used in our analysis the region of interest is just above the centre of the image. The brownish video below is that for AIA193, the bright flash is the tell tale sign of a flare. The brightness of the flare causes very bright saturation and blooming above and below the flare region on the CCD detector and caused extended diffraction patterns to spread out across the SDO imager. The next (yellowish) video is for the AIA171, the red/green mottling shows the magnetic field strength and polarity from measured using the HMI.

The steps followed in our analysis procedure are as follows

Step 1.  Image translation

From our video generated using the output from helioviewer we observe the features moving from left to right as the sun rotates. To perform our analysis we translate the images to correct for this rotation so that features of interest are at a fixed location within each image. From the relatively short  helioviewer movie it appears that the rotation is constant for all locations on the solar disc. Since the sun is a gaseous plasma there is a variation of this rotation with the solar latitude. This is given by the relation.

  A= 14.713 ± 0.0491 °/d
B= -2.396 ± 0.188 °/d
C= -1.787 ± 0.253 °/d

phi is the solar latitude and omega is the solar rotation. Since the solar radius is 695.5Mm and the solar disc in our images is 3304 pixels in diameter, each pixel corresponds to 0.421Mm. Using the solar rotation the translation is approximately
\Delta is the resulting shift in pixels, P_s is the pixel sclae (0.421Mm) \omega is the solar ration in radian/s. t_i and t_0 are the times of the observation and the initial observation respectively. Aschwanden 2011 (ref. 6 ) provides details of this correction in their paper. Our implementation of this correction is in the Matlab script file translatelopsequence.m which we applied to the full sequence of images

Step 2. -Region selection and image cropping

Our region of study was NOAA group 12003 located at N06W37, no flares were reported for this group on this day. We set the challenge of seeing if we could detect changes in the group which occurred as a result of the flare event at NOAA group number 11996  located at N14W76. Inspection of group 12003 indicated that it may be sufficiently near and it had a number of features which we sufficiently stable and well defined to measure.  With the script  croploopsequence.m we applied the Matlab imcrop function to the sequence of translated images.

Step 4 - loop examination and comparison

Tracking the changes in a particular coronal loop is actually quite challenging. The eye is good at picking out features as they change within a sequence such as that within a video. However when we study individual frames and wish to pick out and detect the evolution of an individual loop this is challenging. The process of identifying the loops was achieved by altering the image contrast using the matlab histeq function, for each image set. We compared these equalised images with the original cropped image, with a check against the video frames and by marking off points on earlier and later frames in the image sequence. It was necessary to undertake this process iteratively sometimes returning to correct earlier time steps. In this way there is greater confidence that we are tracking a particular loop.   

For our preliminary examination of the loops we prepared the comparepairs.m Matlab script. This was used to experiment with some of the functions within the image processing toolbox and computer vision toolbox. This included tests of the image registration functions, image registration is the process of aligning images of the same scene. We also experimented with some of the edge detection techniques. These tests were undertaken using the processloopsequence.m, detectfeatures.m and houghlines1.m scripts. Given the dynamic nature of our coronal loop images it was quickly realised that automating a procedure for isolating and measuring specific loops according to a prescription is challenging. It was decided to manually measure a subset of the images. We selected 13 images centered around the flare event. 

For the AIA193 data set the event occured at 22:40:43 on 20140312, this corresponded to image 89 in the overall data set. Our subset was the following images
The fifth image in this subset corresponded to the flare event.
These images correspond to the following time in minutes
20.3   50.7   80.7   95.3  110.7  125.8  140.1  154.7 171.0  198.7  214.5  231.1  261.3

For the AIA171/HMI data set the image nearest to the event was taken at 22:44:23 on 20140312, this corresponded to image 54 in the overall data set. Our subset was the following images
The fifth image in this subset corresponded to the flare event.
These images correspond to the following time in minutes
24.4   54.0   84.4   99.4  114.4  129.4  144.4  159.4 174.4  202.4  217.4  232.4  262.2

For the indices above note that we did not select all the images for the AIA193 we selected those which were closest in time stamp to the AIA171 images

step 5 - loop measurements

For this study we measured the length, height and width of a well defined loop within the group and which remained sufficiently define over our selected  sequence of 13 images. We also measured the loop intensity. The measured parameters are illustrated by the diagrams a,b and c shown below.

We used the  matlab histeq function and balanced the image contrast, it was then easier to identify and measure loop features. The images show curves plotted on the selected loop elements after the initial loop observations, these were the lines which were used to measure the loop intensities. The most effective method is to click and select points from features, for example we clicked points along the flux loop. For this step we used the plotimpair.m script, the measurements were  stored in the loopproperties_img2.m, for use by the later image processing routines. In order to ensure we were measuring the same and correct loop this process was repeated a number of times.

For the image below we have selected the image nearest to the flare event. For the AIA193 we can see clearly the flare event. The brightness of the flare causes saturation and blooming above and below the flare region on the CCD detector it also caused extended diffraction patterns to spread out across the SDO imager. The AIA171 iamge also includes the mottled patches from the HMI corresponding to the strength and polarity of the magnetic field.

The results from this measurement process are shown below. In the following images its important to note that the flare event occured at 110.7mins (for AIA193) and 114.4  (for AIA171)
Variation of Loop Height with Time (Flare event at 110.7 minutes)

Variation of Loop Length with Time (Flare event at 110.7 minutes)

Variation of Loop Width with Time (Flare event at 110.7 minutes)

In the above images it is encouraging that the measurements agree to within the errors we expect for this measurement process. There is a couple of points which are questionable indicating the possibility of incorrect lop/feature selection.

Step 6 Loop thermal analysis

The ability to measure and study thermal characteristics and the thermodynamics of features in the solar atmosphere is a highly valuable tool, this  has developed increasingly since the availability of multispectral, high resolution and high cadence solar imagery. In our final analysis step we used the technique of isothermal analysis to determine the temperature of the coronal loops. This analysis of the loops was computed using the matlab script, loopintensitycomputation.m. In the previous section we had attempted to identify and highlight the full profile of a particular coronal loop. Very often a section of the loop would overlap or intersect with a different loop. For the thermal analysis we used only well defined sections of the loops for the analysis. We edited and stored these sections in the script file loopproperties_partialloops_img2.m. Using these sections of loop we compute the pixel profiles for the gray scale image generated using the the Matlab rgb2gray function. The pixel profile was computed using the improfile function. Using our loop sections as input for the  improfile function we can use the resulting image profile to compute a pixel data number (DN), representing the intensity of the emission from the loop. We describe the procedure for computing the DN values below.

The result of this analysis lead to the following image, here it is important to note that the flare event occured at 110.7mins (for AIA193) and 114.4  (for AIA171)
Measured Loop Intensity and its Variation with Time for the AIA171 and AIA193 Images
In the above image we see that the loop intensity falls initially but that there is rise in intensity occurring after the flare event. To determine if this is at all significant we see if these values give temperature values which are consistent with typical solar coronal loops. To make sense of this it is necessary to recognise that the SDO AIA cameras uses different filters which have varying and known sensitives to radiative emissions of different wavelengths. AIA has seven EUV and three UV–visible channels. Six of the EUV channels are centered on a different spectral line of ionized iron and allow the construction of relatively narrowband temperature maps of the solar corona from 0.4 to 20 million degrees. The plot below (see ref. 12) shows the response curves for the different filters on AIA, these curves show how the image intensity varies with temperature. It can be seen how the different filters have different sensitivities for different temperatures. Using this information it is possible to use our measured intensity data to estimate temperature values.

Response Curves for the SDO AIA Cameras, Responses are Shown for the Different Filters (See ref. 12)
In the figure above the fluxes are given in units of data numbers per second (DN s−1) these counts  are normalised by detector the exposure time, which is different in each filter. Since different spectral lines have varying sensitivies for different temperatures this fact can be used to identify characteristics for different regions (i.e. different heights)  within the solar atmosphere (e.g. see solar facts1). The temperatures for the main spectral lines are given in the table below, the full information of these spectral lines for astrophysical plasmas is provided by the Chianti database

Here is the technique for computing the DN's from ref 12. For each of the loop sections we averaged the loop- and background-pixel values and calculated standard deviations.  A background region is adjacent to the loop under consideration. We subtracted the averaged background results froom the averaged pixel counts for each loop section. Averages were taken over the loop length or the background section. The background-subtracted averages and uncertainties were then normalized by the exposure time for the appropriate filter, this results in units of Data Numbers s−1. Fluxes are now given in units of data numbers per second (DN s−1) after normalizing the counts (DN) by the exposure time, which is different in each filter. Aschwanden (ref 6.)  uses spline fitting for flux profile averaging. For each selected loop they mark the start, midpoint, and end position of the loop spine and interpolate a 2D-spline function with steps of one pixel in order to obtain a smooth curvature. Our implementation is presented in the script file loopintensitycomputation.m the AIA filter characteristics  are stored in loopproperties_partialloops_img2.m. For an isothermal analysis of flux loops we can use the following relation

Flux_171 and Flux_193 are the measured fluxes for the 171 and 193 filter respectively. Resp_171 and Resp_193 are the detector responses for a given temperature for the 171 and 193 filter respectively. Armed with this we can use this ratio to estimate the isothermal loop temperature. Using the detector response values the response ratios for each temperature are shown in the plot below.
Isothermal analysis for the 171-to-193 ratio.
Using the ratio of the fluxes for each of our loops we can read the loop temperature from the above plot (since we know flux ration is the same as the response ratio). This procedure enabled us to produce the following temperature profile. In the resulting image, below, we note that the flare event occurred at 110.7mins (for AIA193) and 114.4  (for AIA171). There is a hint of a temperature rise following the event.
Our approach is a simplistic preliminary attempt and we refer the reader to the excellent work and reviews in references 7,8,9,12,13,22. In this area there is a controversy surrounding the challenge of computing the cross-field temperature distribution for a loop. Also, are loops isothermal or multithermal? Is the observed loop a single flux tube or a collection of tangled magnetic strands? Here, a loop is a discernable structure in an observation and a strand is a fundamental flux tube where the temperature and density are constant perpendicular to the magnetic field. A loop could be composed of a single strand or many unresolved strands. Resolving this controversy has important implications for the coronal heating problem. Results from isothermal analysis, which is the average of the true plasma distribution weighted by the instrument response functions, appears to be deceptively low. These results have potentially serious implications and show substantially smaller temperature gradients than predicted by standard models for loops in hydrodynamic equilibrium these have been used as strong evidence in support of footpoint heating models.

Using the technique of differential emission analysis  has shown that the loop plasma is indeed multithermal and not isothermal. The Differential Emission Measure (DEM) ξ(Te) is a measure of the amount of emitting plasma along the Line Of Sight (LOS) as a function of the electron temperature Te. The intrinsic difficulties involved in this inverse problem lead to many complications in its inference, making its interpretation ambiguous (see reference 12 and 13).


We have demonstrated the application of the matlab image processing toolbox for studying imagery of coronal loops produce using SDO AIA camera.

We studied a loop system in a region close to an M2,5 class flare. Results appear to indicate an increase in the loop intensity after the flare event. This may not be surprising given the  amount of energy liberated by the flare event a and the possibility of radiative heating of our coronal loops.

A study of a single group is a provider of an indication, further results would be required to clarify this conclusion. Measurements of loop systems well isolated from flare events would be required, also necessary is a measurement of loops near and at varying distances from flare events.

It is intriguing that we observe the initial intensity peak for the 171 and 193 these peaks are separated by almost 49 minutes. It is a risk that we have selected a different loop although positional data and the iterative process indicates we've reduced this risk. 

Although we have established a procedure using an iterative approach to make these measurements it is a time consuming process. It is difficult to apply the Matlab visualisation toolbox methods without significant effort and development. A more successful approach might be that of a citizen science approach for example see the zooniverse sunspot project or the zooniverse solarstorm watch. This method beneficial when computer visualisation algorithms find it difficult to pick out and match the details in the images we have studied. Approaches used could be machine learning such as a neural network presented with a training set of coronal loop imagery, the training set includes a sequence of cropped  (and suitably translated images) along with a desired output.

We refer readers to reference 22 for a detailed review of the physics of coronal loops.
Key issues for the future remain the loop fine structure and dynamics. We need to address what are the ultimate elementary loop components, and whether they are unique or determined by local conditions, i.e., what determines the section scale size. We need the highest possible spatial resolution, probably in different bands. Further investigation of the plasma fine thermal structure and dynamics requires also high spectral resolution. High resolution broad-band X-ray spectroscopy is foreseeable to probe the hot components, signature of impulsive heating. Also the investigation of temporal variations still deserves attention. If loops are really so dynamic and subject to a distribution of heating events, whatever it is, the light curves are very difficult to interpret and signatures of any possible small-scale events are confused in a storming activity and by the plasma inertia.


We have used the matlab image processing toolbox to process a sequence of images of solar coronal loops taken by SDO. In hindsight it would have been better to use the VSO, the heliophysics integrated laboratory or myexperiment (ref. 15 and 16) the sunpy or or the SSW. Using this we can extract the image metadata. The matlab toolbox also provides FITS toolbox, this provides highlevel functions such as fitsread and fitsdisp which can be used to read a FITS data file and the metadata respectively.

  1. Getting Started with Image Processing Toolbox-2014b
  2. Image processing examples - matlab central file exchange
  3. Image Processing with MATLAB-Pascal Getreuer
  4. MATLAB 6.5 Image Processing Toolbox Tutorial
  5. Automated Temperature and Emission Measure Analysis of Coronal Loops and Active Regions Observed with AIA/SDO Markus J. Aschwanden, Solar Physics 2011
  6. Solar Corona Loop Studies with AIA: I. Cross-Sectional Temperature Structure Markus J. Aschwanden, Paul Boerner, > astro-ph > arXiv:1103.0228, (ADS 2011ApJ...732...81A)
  7. Coronal Loop Observations and Diagnostics
  8. Evidence for Nonuniform Heating of Coronal Loops Inferred from Multithread Modeling of TRACE Data Aschwanden et. al. 2000ApJ...541.1059A
  9. Time Variability of the ``Quiet'' Sun Observed with TRACE. II. Physical Parameters, Temperature Evolution, and Energetics of Extreme-Ultraviolet Nanoflares Aschwanden et. al.2000ApJ...535.1047A
  10. Solar corona images observed using solar eclipses
  11. A Multi-Wavelength View on Coronal Rain Muller 2004 (see also ADS Muller  20052005ESASP.600E..30M  )
  12. Isothermal and Multithermal Analysis of Coronal Loops Observed with AIA, Schmelz, 2011ApJ...731...49S
  13. On the Accuracy of the Differential Emission Measure Diagnostics of Solar Plasmas. Application to AIA/SDO. Part II: Multithermal plasmas    C. Guennou et al 2012 arXiv:1210.2302 [astro-ph.SR]
  14. JSOC SDO AIA Metdata keywords
  15. Workflows for Heliophysics, Journal of Grid Computing, Anja Le Blanc et al, September 2013, Volume 11, Issue 3, pp 481-503
  16.   PIERANTONI, G. et al. A WORKFLOW-ORIENTED APPROACH TO PROPAGATION MODELS IN HELIOPHYSICS. Computer Science, North America, 15, jul. 2014. Available at: <>
  17. High Speed Coronal Rain Muller et al. 2005A&A...436.1067M
  18. E. O'Shea, D. Banerjee, and J. G. Doyle,  Plasma condensation in coronal loops,  A&A 475, L25-L28 (2007) , DOI: 10.1051/0004-6361:20078617
  19. Ionson, J. A., Resonant absorption of Alfvenic surface waves and the heating of solar coronal loops, 1978ApJ...226..650I
  20. James A. Klimchuk et al., Can Thermal Nonequilibrium Explain Coronal Loops? 2010 ApJ 714 1239
  21. Schrijver, CarolusJ., Catastrophic cooling and high-speed downflow in quiescent solar coronal loops observed with TRACE, Solar Physics, February 2001, Volume 198, Issue 2, pp 325-345, DOI 10.1023/A:1005211925515
  22. Fabio Reale, Coronal Loops: Observations and Modeling of Confined Plasma, Living Rev. Solar Phys. 11 (2014), 4.  doi: 10.12942/lrsp-2014-4

Wednesday, 4 February 2015

Stable Numerical Solution of PDEs Using Hyperdiffusion

The hyperdiffusion method is a technique for the numerical solution of partial differential equations (PDEs) in a way which is numerically stable. It is known that physically interesting problems involve shocked and unstable systems, obtaining stable solutions for such systems may be numerically challenging. We review numerical methods for solving PDEs and  attempt to understand how numerical instabilities arise. This is aided with the presentation of an example program with matlab and scilab to understand numerical instability and the different stabilisation methods.

The method known as finite differencing is used to solve PDEs when analytical methods are not available.The approach adopted is to map out a problem onto a mesh of finite elements. For example, consider the first order differential equation
A number of finite differencing techniques are available for solving this.
The first order Euler method solves this as
y_(n+1) and y_(n) are values of y at the n th and (n+1)th mesh point, h is the difference between the two points

f(x_n, y_n) is the function evaluated at the nth point. To second order the Euler method is written
Here, we have used the notation defined above, a third order Euler method also available. One of the most powerful soultion techniques for first order PDEs is the fourth order Runge-Kutta method, using this method the solution may be obtained using the following expressions.
Now consider more general PDEs of the form
Equation 10: General PDE
Consider a Taylor expansion around the point x, so we may write
Manipulating the first expansion we obtain the derivative of F in terms of the forward difference
the derivative can also be obtained using the backward difference

Taking the difference between the two expansions, gives the derivative in terms of the central difference
Notice that we have neglected the higher order terms in the Taylor expansion. Central differencing formula for the derivatives to order h^2 are as follows:

The notation used above is as follows
and for example
The central differencing expressions to order h^4 for the first three derivatives is as follows

Now, solve equation 10 using first the forward difference
and then a backward difference
Combining the forward and backward difference we can solve the PDE using the following
This method is known as the Crank-Nicolson method. It is a so called implicit method in that knowledge of the solution at time t_{i} gives a relationship with the function at time t_{i+1}, which is not given explicitly. We now apply these methods to the solution of the following advection problem.
Equation 29: Advection problem
For the above equation, when c=u, then we have an equation called the Burgers equation, this so called inviscid equation can generate solutions which can develop discontinuities and shock waves. In our matlab and scilab script we attempt to solve this numerically for a propagating square wave and a sine wave. Solutions for the central differencing method at two different times are shown below.
Advection of a square wave at time t1
Advection of a square wave at time t2
The first image shows the initial square wave and how it is advected after a small number of time steps, it can be seen clearly that numerical disturbances arise at the leading and trailing edges of the wave form. At the even later time t2, the noise has grown and it is apparent that these solutions are numerically unstable. Why is this?

To obtain stabilised solutions a number of solution approaches are possible including the hyperdiffusion approach, but first we consider another approach to obtaining numerically stable solutions by considering the Lax-Wendroff method. Consider the following PDE.
To obtain numerically stable solutions of such a problem we may use two stages. In the first stage the solution is advanced by half a time step and by half a mesh point in the forward/backward direction as follows

We advance the solution using the results from the first step as follows
Using this approach our advection problem, equation 29 may be solved using the following finite difference expression.
Using VonNeumann stability analysis the stability of the solutions can be checked, with solutions of the form
For the lax-Wendroff method the amplification factor A is found to be 1. If A<1 the amplitude for numerical diffusion is found to decrease. If A>1 the amplitude for numerical diffusion is found to increase and the solution is numerically unstable. The result of the application of this solution method is shown below.
Advection of a square wave at time t1 using the Lax-Wendroff method

Advection of a square wave at time t2 using the Lax-Wendroff method
Using the matlab script file comment out the lines for the central differencing method and the hyperdiffusion corrections, Also comment out the correct lines before the statement, for i = starti:finishi. It is necessary to set the correct values for starti and finishi, for the Lax-Wendroff method we use the second order differencing values. When this script is run we should observe the advection of the square wave as illustrated by the two figures shown above. Once again, observe the numerical instability at the leading/trailing edges of the square wave. However at the later time t2 although there is still numerical noise this has not grown significantly and we can successfully compute solutions at later times.

We now consider the source of the numerical noise and how the hyperdiffusion method may be used to control that noise.

To help understand the source of the numerical noise consider the solution of the following PDE
let u be cosntant and solve the problem in 1D, so that we write 
Using a Taylor expansion we can write the follwoing Taylor expansions
Taking the difference of these two expressions our finite difference version of the PDE becomes
If we choose a solution of the form
then we can obtain an expression for the frequency as follows
It can be seen that the second term in the series causes a reduction of the frequency and the solution will change shape. The phase speed of the signal is
Shock waves and square waves may become dominated by high frequencies. Therefore, f we can use a method, such as central differencing to truncate the solution and apply this near the discontinuity it may turn out that the solution is actually dominated by the higher frequencies. An nth order term may be written as
For the 1D problem the second order correction may be written as
The introduction of the terms such as this leads to a damping term which can be used to damp out some of the higher frequency terms of a given wavelength. Solving the above equation solutions of the following form are obtained
The coefficient nu_2 damps out the spatial wavelengths.

Now consider the 1D problem
In this problem we have numerical discontinuities arising from shock, advection and sound propagation. If the introduced term is adjusted to balance the jumps in the solution arising from the higher order terms then a solution may be stabilised. Following the definition of our damping term earlier a general expression for the viscosity is as follows
The problem now is to set the coefficients so that the discontinuities are successfully cancelled. So consider our solution of the form
If a wave of length 3delta x dies out quickly then after 4 lengths it should be cancelled.
The time t_e is the time such that the signal will decay  to 1/e of its original what must nu_2 be in order to satisfy this.
If we define
Now, define
Using these definitions we obtain

Thus our differential equation with numerical correction becomes
We have computed the value of v_2 to enable the correct cancellation of the higher order terms. We can compute the wave speed and we can compute delta x from the Courant condition. Now consider the equation of momentum conservation from our perturbed system of equations for ideal MHD. Below rho_b is the fixed background density and rho is the perturbed density. For full details see reference 1 below.

We will now consider how we compute the hyperdiffusion corrections for this system of equations
The above expression takes the same form as the corrective term we described earlier, the second term is a correction for the shock viscosity.

delta^{3} is a third order difference and delta^{1} is a first order difference, i.e. using central differencing,

These expressions will be used to pick out the numerically unstable parts of the solution for our computational mesh. For our MHD code we split these contributions into a left handed and right handed differencing contributions.
We express this numerically as
The maximum is taken over the left and right differences over the computational mesh. We have implemented this term in the example matlab and scilab codes, here each of the terms are named d3l,d1l,, d3r and d1r.

The shock viscosity term is calculated for negative compression, the compression term is 

Only regions undergoing compression are corrected using the shock viscosity. The hyperdiffusive operator 
c_s is the sound speed, v_a is the Alfven speed.
is the hyperdiffusion function and

We may write
These terms have been implemented in our matlab/scilab code through the terms labelled nushk(i) for each of the mesh points.
Advection of a square wave at time t1 using the hyperdiffusion method
Advection of a square wave at time t2 using the hyperdiffusion method
The pictures above illustrate the advection of a square wave using our matlab/scilab code with the implementation of the hyperdiffusion method described here. We can see that the correction improves on the Lax-Wendroff approach and is effective at preserving the wave shape.

We now consider the hyperdiffusion terms for the equation of continuity and the evolution of energy.
For the energy term

Where the thermal energy perturbation is used and given by,
In the expression above, epsilon_klm is the levi-civita symbol