MA2895 Assignment: Corrections, Comments and FAQ

Links to the FAQ.

Access to h: LaTeX Finding roots of cubics tasks Partly graphics task FDM task

Links to the comments.

General comments cube_root_fun cardano_fun test_cubic_iter cubic_plot basic_fdm numerov_fdm test_fdm

Corrections

  • In the 17th Feb 2020 version the table briefly describing each of the files the m-file in the graphics task has a name ending in _cubic_plot.m and not _gr_task.m (as given in the table). The file that you were sent has a name ending in _cubic_plot.m and the _cr_zip_file.m file expects this. You were sent the corrected version of the entire pdf on 20th Feb 2020. This remains the only correction to report.

    Answers to FAQ

    Access to my h: drive files using VPN

    The starting point to set this up is to follow the link below. You will need to login with your usual Brunel details.

    Click here for the Brunel connect portal http://connect.brunel.ac.uk

    For a VPN connection you need to click on the box to the left of ``AnyConnect'' and then click on Next. If you are doing this for the first time then you will need to click on ``New Device'' to get a device in this list. I am not in a position to test this at the moment but from what I can remember it is at this stage that you will download and install the software on your computer.

    After you have registered your device this should appear in the the list of devices.

    Please always check what is currently said although in case it helps the link (click here) contains what it showed for me on Fri 27th March 2020.

    From the window which had the ``AnyConnect'' link mentioned earlier there is a link with which says ``Brunel Files''. You should click on that and follow the instructions to access your h: drive. The instructions here depends on the user.

    If you need further help then please send an email to 

     CEDPS-IT@brunel.ac.uk

    The file print2pdf.m and other files in h:\mybin20

    The file print2pdf.m is located in 
    h:\mybin20
    if you did the initial set-up instructions. If you cannot access the h: drive now, or it is more difficult, then a zip file called mybin20.zip containing most the files is available if you     click here.

    The finding roots of cubics task

  • What should he shown when the input argument called out is different from 0? Also, does it matter if my numbers are not identical with those in the assignment document?

    I do not specify exactly what should be shown but ir should give an indication that the iteration is converging. In $b$ is given and the function is $\phi(a)=a^3-b$ then showing that the norm of the $\phi$ values is decreasing rapidly is sufficient. In the test with 101 values on a circle I do not want to see 101 values for each stage of the iteration. It is your choice whether or not you should the values when the number of values is small (say 1 or 2).

    The letters used in the description vary from student to student. If I use a_0, a_1, ... here then a_0 corresponds to b^{21/64}, a_1 is from the Newton step and a_2, a_3, ... are from the even faster method.

    When numbers are small you should use the %e format. If you use the %f format then it will probably just show 0 which is not a satisfactory display.

  • It does not matter if the 3 roots generated by the cardano function are in a different order to what I show. The first root should however be real.

  • If w is the vector with the 3 roots using the cardano method and w_e is the vector with the 3 roots using the Matlab function roots() then one of the tests is to check that they are the same 3 values (approximately). It is unlikely to be sufficient to just have w-w_e and also you need to take care that both w and w_e have the same shape. You can ensure that they have the same shape by having w(:)-w_e(:) but you will only get close to the zero vector if the roots from both methods are in the same order. For each w(k), k=1,2,3 you need to find the nearest entry in w_e when you compare.

    The partly graphics task

  • The formula for the 2 areas is the unnumbered displayed equation before equation (11).

  • There is information about graphics commands in section 7 of the Matlab revision notes. The only Matlab commands that I have in my solution of each version of the assignment are the following. 
    axis equal
clf
figure
fill
hold off
hold on
plot
title
    I do also use the following when I am running many scripts. 
    close all
print
    the command ``close all'' closes all figures and ``print'' is used to get a hard copy version. I used print to generate the .eps file that you were sent and used in the file print2pdf.m to help generate the .pdf graphics file that you were sent.

    All the above command are mentioned in the revision notes.

    The finite difference methods task

  • The main difference between the numerov_fdm function and the basic_fdm function is that the tri-diagonal matrices are different and the rhs vectors are different. With the Numerov method the functions q and r are evaluated at all the points in the vector x whilst in the basic method there is no evaluation at the points a and b where [a, b] is the interval. In the first argument to spdiags the quantity has 3 columns with the first column containing entries on the subdiagonal and the third column containing entries on the superdiagonal. In each of these cases one of the entries is not used and you need to know which one to be sure that the others are in the correct position. There is some information about this in the Matlab revision notes and there is a section in chapter 4 of the MA2715 lecture notes which deals with this.
  • When the program breaks down at the statement involving backslash to solve the linear system this is often because the rhs vector does not have the correct shape, i.e. it is not a column vector with the same number of rows as the matrix. The correction needed here may not be in the function file which has the backslash statement. The rhs vector depends on the function r(x) and this in turn depends on uex and uexdd both being done correctly in a vectorised way. As a first check you might try the following statements. 
    x=0:0.25:1
x=x(:)
uex(x)
for k=1:5
  uex(x(k))
end
    Do you get the same 5 entries as column vector two-times. Then try 
    r(x)
    Does this also give a column vector of length 5.

    One way to check if uexdd is consistent with uex is to have the following. 

    x=linspace(a, b, 11);
h=eps^(1/3);
uddapprox=(uex(x+h)+uex(x-h)-2*uex(x))/(h*h);
    You then compare uddapprox with uexdd(x). The finite difference approximation will be close when uexdd is correct. If norm(uddapprox-uexdd(x), inf) is large then it likely that uex and uexdd are not consistent.

    General comments

  • The following is a general comment about doing division correctly in a vectorised way. Suppose you wish to evaluate $\sin(x)/x$ at a few points. The following script indicates how to do this and some similar attempts which give a different outcome. 
    format compact
x=[0.1, 0.2, 0.3]

y1=sin(x)./x

y2=sin(x)/x

y3=sin(x)\x
    The output generated is as follows. 
    x =
   0.100000000000000   0.200000000000000   0.300000000000000
y1 =
   0.998334166468282   0.993346653975306   0.985067355537799
y2 =
   0.988380498729264
y3 =
                   0                   0                   0
                   0                   0                   0
   0.338386336182412   0.676772672364825   1.015159008547237
    The output given by y1 is what is wanted. More complicated computations are taking place in the case of y2 and y3. If you type 
    doc slash
    then this will explain what happens in the y2 and y3 cases.
  • Take care to distinguish between the digit 1 and the lower case letter l. 1e-6 is for 10^{-6}. In Matlab if you have le() then this is a function version of the <= comparison and the function need 2 arguments.

    Comments relating to the cube root function file

  • Take care with the stopping criteria in the iteration to get the cube root. Many of the issues that I have seen relate to a set-up such as the following. 
    format compact

b=[9, -1+1i];

phi =@(a) a.^3-b;

a3=[2, 0.8+0.8*1i];

vrow=phi(a3)

vcol=phi(a3(:))
    When you run this you get the following. 
    vrow =
  -1.00000 + 0.00000i  -0.02400 + 0.02400i
vcol =
   -1.00000 +  0.00000i    9.00000 -  1.00000i
  -10.02400 +  1.02400i   -0.02400 +  0.02400i
    The different output is because the vector b used in the definition of phi has the shape of a row vector. The evaluation of phi needs to be with the vector having the same shape as b in your cube root iteration. Thus in the above the first version using the row vector a3 is likely to be what is intended. With releases of Matlab within about the last 2 of 3 years the other version with a3(:) no longer gives an error but creates a matrix, i.e. it does not breakdown with a column vector minus a row vector. The norm of vrow and vcol are thus different.

    To deal with any shape of b, and assuming that a3 has the same shape as b, to get a column of entries for the function evaluation it is best to put 

    vcol=phi(a3);
vcol=vcol(:);
  • Take care to distinguish between the maths description and what you put in Matlab statements. In maths we write statements such as a_{k+1}=a_k-\phi(a_k)/\phi'(a_k). This is to distinguish between a_k and the next value a_{k+1}; In a program you might do this with statements of the form 
    v=phi(a);
dv=dphi(a);
a=a-v./dv;
    The above assumes that functions phi and dphi have already been set-up. Note also, you do not have to put everything in one statement. In the version above the term v can be used in a test without needing to use phi() again.
  • In the starting version of the file that you are given there are statements to create a rescaled version of the number whose principal cubic root is to be determined. The iteration part works with the rescaled version. The statement to finally set the value to return involves multiplying by 2^.e and this needs to be at the end of the function file.

    Comments on the Cardano function file

  • In part of the description it is stated that you need to set something with two possibilities for the value and the decision is given as to which you should pick. I have seen more than one person do this incorrectly. If the 2 possibilities are denoted by w1 and w2 then you need of one of these. Some people are setting the value to one of abs(w1) and abs(w2) which will be correct sometimes and not correct in cases when the correct answer should be negative. To state this another way, the set containing w1 and w2 is not always the same as the set containing abs(w1) and abs(w2).
  • The first root that you create in the cardano function must be real if you are following the algorithm. If you find that when you are testing the function that the first root is not real in some cases then you may still have gathered the 3 roots but you have not followed the algorithm specified.

    Comments on testing the cube_root and cardano script

  • The only loop you are likely to add to the file beyond what was in the starting version is a loop to consider each of the 3 radius values.
  • With a vector t already given with 101 values you can get 101 values on the unit circle with exp(1i*t) or ( cos(t)+1i*sin(t) ). Here 1i is the pure imaginary number which is the square root of -1.

    Comments on the cubic plot script file

  • If your plot in figure(id3) looks different to the one in your assignment document, and which was also sent to you as a separate file to be able to use in your ma2812 latex report, then this suggests that the constant term d in your cubic is not set correctly. The statement for d to use is given in the first item of section 4.2. This depends on dfac, which was already in the file, and it depends on beta and gamma being done correctly. Remember that this is an individual assignment and the precise description and the precise instructions vary and the variations were randomly generated.

    For a plot of a circle to look a circle you need the statement 

    axis equal
  • Most of the things to do in this task are connected with the plot associated with figure(id3). Immediately after the figure statement I recommend that you put clf. As I mention, you will want to use the ``hold on'' and ``hold off'' mechanism if you are attempting to do all parts. The ``hold off'' statement will need to be after the last plotting statement that you use. If this is not done then the previous plotting parts may appear very briefly but then are replaced by the plotting instructions put after ``hold off'.
  • In the task to generate the smaller circle I put (with latex syntax) x_m to refer to the middle of two of the roots. In your program you might put xm or x_m for a corresponding variable name. There is no reason to create a vector called x and store the mid-point in one of the entries. Similarly when I put (with latex syntax) y_c the notation is just indicating the vertical component of the centre of the circle (with abs(y_c) being the radius in the context of the task) you should not be creating an array y and referring to the entry in position c. You can use yc or y_c or any other sensible name for a variable for this.
  • The maximum number of possible terms in the loop in the algorithm to get the smaller circle is not specified. A number such as 10 might be used. If everything is done correctly then the loop should finish well before 10. The first value of yc to try will be too big but visually it would likely that a circle of 1/2 or 1/4 or 1/8 of this should fit into the shaded region.

    Comments on the basic fdm function

  • You will not know if the function works until you have statements which use it and such statements can be part of the test_fdm script file.

    If you are not ready to do the part involving creating uexdd for your specific problem then you can use other simpler functions.

    Both methods are exact when uex is just x^2 and hence you could use the following in a test. 

    a=-1;
b=2;
q =@(x) ones( size(x) );
uex =@(x) x.^2;
uexdd =@(x) 2*ones( size(x) );
    The Numerov method is also exact when uex is x^4 and thus you can check with the last 2 lines replaced by the following. 
    uex =@(x) x.^4;
uexdd =@(x) 12*x.^2;
  • When c is a column vector and A is a square matrix with the same number of rows then you solve a system with 
    A\c
    Note that it is backslash.

    Comments on the Numerov fdm function

  • The statement above about testing basic_fdm applies also to testing the numerov_fdm function.

    Comments on the test_fdm script

  • In some attempts so far there have been mistakes in the expression for the 2nd derivative whereas everything else has been okay. It is not difficult to get the second derivative. If you use any software to do this for you then you will not get the correct expression if that software has done some rounding operations to neatly present the answer.