These are the first INSTRUCTIONS to my MatLab demo programs.
The language is
ENGLISH because of my cooperation with American colleagues.
I regret that some interface-commands do not work in the UNIX-version
as well
as on my MAC.
My demos are written with a double purpose; I hope they are
mathematically
instructive and thus support our first year courses; I also hope
that they help
to get started quickly with programming in MathLab.
*** The IDEA is to run the demos AND to look at the code AND to
use MatLab help
(for example type: helpwin in the command window to get a
browser into help)
AND to read the following comments.***
If you give in UNIX the commands: matlab_getdemo and then:
matlab
then a matlab folder is created in your directory, my demos are
copied into it
and the matlab search path is set to include these. If you want
to modify these
programs you have to use an editor, save your versions under a
different name!
A few general remarks before the first example:
BASIC FACTS, *** READ (0)-(6) before trying anything.
***
(0) MatLab has a numerically based portion and a symbolically
based part. My
demos illustrate the numerical part.
(1) The programs which one writes in MatLab are saved with extension
.m and therefore
called m-files. They come as "script" and as "function".
Every variable which gets
a value assigned in a "script" still has its last value
when running the script has
ended. Every variable, which is used inside a "function"
and does NOT occur in the
input or output list, no longer exists after the "function"
has returned its output.
Since MatLab-functions can return lists of output, they replace
procedures,
subroutines and functions (of other programming languages).
(2) MatLab uses automatically COMPLEX NUMBERS, but recognizes
real numbers
so that zero imaginary parts are usually not printed; see Nr3Complex
(3) MatLab has a lot of LINEAR ALGEBRA built in; this is built
on some efficient
abbreviations. It is not useful to read MatLab programs without
knowing about
these abbreviations. (More in Nr2Maps2D). TRY the following NOW
and type (in the
COMMAND window):
t=[0:10] RETURN-key
This puts the row of integers 0, 1, ..., 10 into the variable
t without the need
of any declarations. This list can be scaled, shifted and put
into functions, TRY:
x=t*pi/10, y=1+sin(x), or x=t*pi/10; y=1+sin(x); <==
Note the ; ;;;;;;;;;;;
In the first version all the numbers are printed in the command
window. In the
second version, with semicolons(;) this output is suppressed!
;;;;;;;
r=[t,t] (or [t t]) makes a row twice as long, c=[t;t] makes two
rows. Note that
the []-brackets are only needed as input help, [1,2,3] and [1,[2,3]]
create the
SAME input, while in mathematical set theoretic notions M={1,2,3}
has three elements,
but M={1,{2,3}} has only two elements, 1 and {2,3}. (This agrees
with the later use
of { } in MatLab.)
cc=c' makes of the two rows of the above c the two columns
of cc.
Try simple examples now, copy each line into the COMMAND WINDOW:
m=[1 2 3; 4 5 6], m'
a=[1 2 3.5; 1.1 2.1 4.1; 1.3 2.3 4.3], b=[0 0 1; 3 7 1; 4
6 5]
plot(a, b), axis([0,5,-1,8])
plot(a, b,'o'), axis([0,5,-1,8])
help plot
M=[[m,[7,8]'];[9 10 11 12]]
More of this in the beginning of Nr2ReadMe.
(4) The graphic default capabilities of MatLab are convenient:
t=[0:10]; x=t*pi/10; y=sin(x); plot(x,y)
shows a graph of sin on the interval [0, pi]. Then type: plot(x,y,y,x)
it shows two graphs, namely of sin and arcsin.
(5) HELP should be used routinely when one looks at demo programs.
TRY: help plot .
HELP can be used on your own programs, for example on my demos.
Help name gives
the first % comment lines of the program (script or function)
called: name.
Initially, I found the ERROR-messages mostly helpful; however,
a missing closing
bracket ] may not be detected and a strange message points to
the next line.
(6) The simplest command to make an INTERACTIVE program is PAUSE,
a big P appears
on the screen and the computers waits for a key of the keyboard
to be pressed.
*Unfortunately, this does NOT work in the UNIX-version*, the m-file
unixpause.m
prints instructions in the command window.
FirstConventions is a program which suggests more such experiments.
It should have
many pause-commands, but since these do not work in the UNIX
version, the whole
program runs once. The arrow-keys then allow to repeat any
of the commands. Since
FirstConventions is a script, one can use all previous variables
in further
experiments.
** The simplest demos are in Nr1RealApp. **
*TYPE the name* of a demo in the command window, RETURN-key to
run it.
Use help if you do not know the number of arguments. Use ;
to suppress output.
tangdemo(2) shows a plot of sin and a tangent on a smaller interval
around 2.
tayldemo(2) does the same with the second taylor approximation;
it illustrates
also the use of FIND and FLOOR and shows how to change colours
in PLOT.
Observe that t2=t.*t squares EACH entry of t and stores it
in t2. Try:
t=[-5:10]/5; s2=t*t', ss=sum(t.*t), plot(t,t.*t)
The output s2, ss is written into the command window; the output
t is suppressed.
taylorap calls tangdemo(a) and tayldemo(a) inside a WHILE-loop.
This illustrates
the graphic input command GINPUT . The computer waits, showing
(unfortunately
only after touching the mouse) a movable cross in the graphics
window.
Press any keyboardkey, then the coordinates of the cross and
the ASCI-code
of the pressed key are returned to the program. This WHILE-loop
only terminates
when the RETURN-key is pressed.
W a r n i n g: The RETURN-key returns nothing from GINPUT,
technically the
empty set [], therefore the question ISEMPTY has to be used
in this program.
Try (by copying to Command Window): a=1, b=[], isempty(a),
isempty(b)
y=sin35(x) is an approximation of sin which combines a taylorapproximation
with
the formula sin(3x)=sin(x)*(3-4*sin(x)^2) (a^2 means a*a).
Note that polynomials can be defined by giving their coefficients,
p=[1,2,3] denotes x^2 + 2*x + 3. Use: help polyval
sinap shows plots of sin, sin35 and (sin - sin35). It illustrates
how to get
TITLES into graphics windows and it shows how to choose AXIS
for PLOT instead
of simply accepting the default values of the axis ranges.
PAUSE (and unixpause) is used.
y=ex43(x) approximates exp by combining a taylorpolynomial with
exp(x)=exp(x/n)^n
Try to modify ex43.
expap compares plots of exp and ex43. The use of STRINGvariables
allows to
choose approximations from below or above. TYPE (to try out
strings):
m=['a', ' b', '..c'], m', length(m), m(3)
s1=['a','b','c']; s2=['d','e','f']; M=[[s1;s2],[' ';' '],['g','h']']
derivdemo compares difference quotients and derivatives. POLYNOM-commands
are used and PAUSE. Change the difference in the quotients!
Why is the symmetric
difference quotient (f(x-h)-f(x+h)/2/h so much closer to the
derivative than the
other difference quotients? Can you estimate
| difference quotient - derivative | ??
Define before running: global h, h=0.3;
then this stepsize h can be changed before each run.
sinap2 shows the approximation of sin by a series in powers of
x*(1-x/pi):
SUM_k {a_k * (x*(1-x/pi))^k}.
The series agrees with more and more derivatives of sin at
the TWO points 0,pi.
The convergence is so rapid that an irrationality proof for
pi can be obtained.
Try to compute a_1,a_2,a_3 so that the above sum has the first
three derivatives
at the TWO points 0,pi the same as sin. Observe the excellent
approximation.
The text written with LEGEND into the window can be moved with
the mouse.
The program writes alternatingly into two graphics windows,
this uses the
command HOLD ON. ( help hold uses the word "toggle",
this means: change the
active state, namely HOLD ON or HOLD OFF, to the other.) PAUSE
is used, but
does not work well under UNIX, I changed it to UNIXPAUSE.
y=pow(x,n) computes integer powers effectively (if ^ is not on
keyboard).
Uses SWITCH, FOR, DEC2BIN
cd=contFrac(num,den,n) computes continued fraction of num/den.
Can either be
called AS: cd=contFrac(pi,32), OR: cd=contFrac(1234567891,
345678321,20)
[x,xf]=contF2dec(cd) is the inverse of the previous. class(x)=double
and
xf is the string numerator/denominator.
differentiable(choice) shows with SUBPLOT several graphs of 1/exp(1/x^2)
and also x^2*sin(1/x^2) on several scales.
w = Interpolate(t,arg,fct); shows the quadratic 3-point interpolation
and
also the fractional linear one.
TRY: Interpolate([0:20]/10, [0.5,1,1.2], '1-cos(x)');
The interpolation subfunctions are also used to interpolate
the inverse
function. Note that functions can be given as input.
w = demoInterpol(t,arg,fct,gct); shows how the interpolation subfunctions
can be used to approximate *pairs* of functions, i.e. parametrized
curves.
TRY: demoInterpol([0:20]/10, [0.5,1,1.5], 'cos', 'sin'); OR:
demoInterpol;
x0 = regulaFalsi(a,b,'fct'); computes a zero in [a, b]. OR: If
'fct' has the
same sign and a, b then a short search for a sign change is
done (this part
mainly deals with input errors). In prinipal this simplest
root finder
replaces the function by secants. Its usual weakness is that
only one endpoint
of the iteration-intervals converges to the zero. This is avoided
fairly well.
TRY: regulaFalsi(0.5,10,'exp(x)-10','graphics'); regulaFalsi(-0.5,10,'sin');
The progress of the iteration can be shown graphically.
t = ratInverse(a,b,w,fct,choice); finds t with fct(t) - w = 0
in [a,b].
This is a strong improvement over the secant method to find
a zero, because
the inverse function is approximated by 3-point-interpolations.
If the points
are far apart, then these interpolations are not monotone.
This has to be
avoided during the iteration; if it occurs, interval-halfing
is used. Once the
interval gets smaller the iteration is very fast since all
3 points converge;
the progress of the iteration is shown graphically,
TRY: t = ratInverse(-1,15,0.5,'sin','graphic');
or t = ratInverse(0.05,20,-0.5,'log(x)');
This program is heavily commented.
pc=cubic(a,b,fa,fb,dfa,dfb) determines the cubic polynomial pc
which has
at a resp. b the values fa, fb and derivatives dfa, dfb.
Illustrates the
multiplication of polynomials with CONV. This is an important
demo; it is
written to accept columns of arguments and then returns several
polynomials.
Simplify the program so that it accepts only numbers (not columns)
as arguments.
Write your own quadratic interpolation pq=quadratic(a,b,fa,fb,dfa).
TRY (by copying into the COMMAND window):
pc = cubic([1;2],[2;3],[1;4],[4;5],[1;2],[2;3]), t=[0:20]/5;
y1=polyval(pc(1,:),t);y2=polyval(pc(2,:),t); plot(t,y1,t,y2)
or try:
pc = cubic([0;1],[1;2.4],sin([0;1]),sin([1;2.4]),cos([0;1]),cos([1;2.4]));
t=[0:30]/10;y1=polyval(pc(1,:),t);y2=polyval(pc(2,:),t); plot(t,sin(t),t,y1,t,y2)
.
democubic shows how well the 11a two-values-with-derivatives approximation
approximates sin on a rather large interval, and arcsin --
AFTER you make the
interval a bit smaller (to get away from the vertical tangent
of arcsin).
Note that with the default [a,b] the interpolation of arcsin
is NOT monotone.
Change democubic so that a,b can be changed with the mouse,
as in 3.)taylorap
by using GINPUT.
spline5 computes a degree 5 polynomial with given values and derivatives
at three points. TRY default: spline5 which approximates sin
OR: see help spline5 to demo the 6 basis polynomials.
splin2ac computes a degree 5 polynomial with given values, 1st
and 2nd
derivatives at two points. TRY: splin2ac.
arg = spline_newton('x^2+1','2*x',1,8,5) finds in the interval
[1,8] the
preimage of 5 for the function 'x^2+1' (also 'x.*x+1'); similarly
arg = spline_newton('exp','exp',-1,5,15) finds the preimage
of 15 under exp.
Programs like this should accept function names as input. The
program
spline_newton shows what commands one has to use so that MatLab
accepts the
two types of function names given as examples. W a r n i n
g: In the DEFINITION
of spline_newton one needs as the last variable varargin which
does NOT(!) appear
when one CALLS this routine. After this preparation the call
of FEVAL evaluates
such functions. -- The iteration produces subintervalls such
that the function
continues to have opposite signs at both ends. Since the double
tangent cubic
approximation of the inverse function is not necessarily monotone
the routine
rejects an "approximation" which is outside the previous
intervall. Also it
avoids that all the improvements happen only to one end of
the intervall.
Check, how much faster than Newton's method this procedure
is by writing your
own Newton method.
The program uses IF ELSE END. Note: if (a==b) needs == ,
NOT = .
Note also: b=(1>0) prints out 1, but b is not a number but
a BOOLEAN variable
which is only TRUE or FALSE; in numerical computations TRUE=1,
FALSE=0.
psq = polyadd(p,q). Note: Polynomials of different degree are
rows of different
length, therefore psq = p+q gives an error.
polytaylor(p) shows different Taylorapproximations of the input-polynomial
p.
(Again: p=[1,-3,7,-5] is interpreted by polynom-commands as
1*x^3 -3*x^2 +7*x -5)
We use again graphic input GINPUT in a WHILE-loop to choose
the position
of the new tangent with a moving cross, then one may press
0 - 5 to choose the
order of the taylorapproximation. Since the order of the approximation
should
be recognizable by the colour we show how to use numbered colours
in PLOT.
Try in the COMMAND WINDOW first p=poly([0,.5,1.5,2.5,3])*10;
this makes
a polynomial with the given ROOTS. Then use polytaylor(p)
and compare p
with several Taylor-approximations at different points (read
title).
[mx,my,r]=contactCircle(a,fa,da,ba) makes [many] osculating circles
of graphs
TRY: x=[-20:20]/10; plot(x,x.*x), contactCircle(1,1,2,2);
contactCircle([1/2,3/4,1],[1/4,9/16,1],[1,3/2,2],[2,2,2]);
simpson_quad and demo_quad are essentially the same numerical
INTEGRATIONS
(called quadratures in MatLab). This is another routine which
should accept
FUNCTION NAMES and indeed does it exactly as in spline_newton.
Recall the
W a r n i n g: The definition needs ..,varargin) the call does
NOT.
TRY as: [mid_quad,sec_quad,simpson] = simpson_quad('x.*x.*x.*x',0,4,4);
[mid_quad,sec_quad,simpson] = demo_quad('x.*x.*x.*x',0,4,4);
simpson_quad gives the three simplest numerical integrations;
in demo_quad
several PLOTs are added to the integration routine, showing
another convenient
default property of PLOT; to understand them one must TRY first:
a=[1 2 3.5; 1.1 2.1 4.1; 1.3 2.3 4.3], b=[0 0 1; 3 7 1;
4 6 5]
plot(a, b), axis([0,5,-1,8])
TRY: help demo_quad. Also, demo_quad calls simpson_quad with
a finer subdivision
and explains how to get from it an even better value by "Richardson"-extrapolation.
order8 = higher_quad(fct,a,b,n) integrates polynomials to 7th
degree correctly
from a to b with n subintervalls. Function-names as in 12.
& 14. accepted.
FOR i=1:n ... END is used.
NOTE: If the support points and the weights of a numerical
integration are
symmetric with respect to the intervall-midpoint and if all
even polynomials
to order 2k are integrated exactly, then automatically all
odd polynomials
to order 2k+1 are integrated exactly. In this order8-quadrature
we chose
four parameters (3 weights, 1 support point position) to integrate
the even
polynomials of order 0,2,4,6 correctly.
cycloids draws, with SUBPLOT, hypo- and epi-cycloids.
TRY: n=5; r=0.4; cycloids(n,r); OR cycloids;