Plots are easy to create in MATLAB. A simple example is a 2D plot of the sine function. We want to plot {{{sin(x}}} for a range of //x// values. First we define the domain by \n\n{{{>>x=linspace(0,10,20);}}}\n\nThis produces a 1x20 set of values from 0 to 10. Note the semicolon will hide the values and not display them. The same set of valuescan be obtained using\n\n{{{>>x=0:1/1.9:10;}}}\n\nTyping\n\n{{{>>y=sin(x);}}}\n\ndefines the vector of 20 values of {{{sin(x)}}} and stores them in {{{y}}}. This is not necessary for plotting.\n\nNow we plot {{{y}}} vs {{{x}}}:\n\n{{{>>plot(x,y)}}}\n\n[img[Sine Function|plotsin.jpg]]\n\nYou can render the image with other attributes. You can make the plot red:\n\n{{{>>plot(x,y,'r')}}}\n\n[img[Sine Function in red|plotsinred.jpg]]\n\nOther colors are possible using the following:\n\n|b|blue|g|green|r|red|k|black|\n|c|cyan|m|magenta|y|yellow| | |\n\nOne can also pick different linestyles:\n\n|:|dotted| - | solid| -. | dashdot| --| dashed|\n\nA plot with the dash-dot linestyle:\n\n[img[Sine Function with dash-dot style|plotsindashdot.jpg]]\n\nOne can plot the values as data points by using various markers:\n\n| . |point| o|circle| x|x-mark|\n|+|plus|* |star|s|square|\n|d|diamond|p|pentagram|v|triangle (down)|\n|^ |triangle (up)| < |triangle (left)||>|triangle (right)|\n|h|hexagram| | | |\n\nThus, a plot with red circles is obtained with\n\n{{{>>plot(x,y,'ro')}}}\n\n[img[Sine Function with markers|plotsinredcirc.jpg]]\n\nOne can also add text to the plots. The common texts needed would be the axis labels and titles. While these can be added to individual figures using the menus on the figure, they can also be put into the script.\n\nFor our plot, we can add\n\n{{{>>title('Plot of sin(x)')}}}\n{{{>>xlabel('x')}}}\n{{{>>xlabel('sin(x)')}}}\n\n[img[Sine Function with labels|plotsinlabels.jpg]]\n\nWe can also plot several functions. Let's plot markers for the above data and plot this with a more detailed version of the sine function with 200 points. This can be done using the folowing set of commands:\n\n{{{>>x=linspace(0,10,20);}}}\n{{{>>y=sin(x);}}}\n{{{>>X=linspace(0,10,200);}}}\n{{{>>Y=sin(X);}}}\n{{{>>plot(x,y,'ro',X,Y,'b')}}}\n{{{>>title('Plot of sin(x)')}}}\n{{{>>xlabel('x')}}}\n{{{>>xlabel('sin(x)')}}}\n\n[img[Sine Function with labels|plot2funcs.jpg]]
This is an introduction to ''MATLAB'' for students. It provides some of the basics needed to get started with ''MATLAB''. This should enable students to be able to write simple scripts and functions for mathematics and physics. Links to other online tutorials are also provided.\n\nThis site was buit with [[TiddlyWiki|http://www.tiddlywiki.com/]].
There are many advanced ideas, which in time might be added to this section.\n\n#Variable Types\n#I/O Commands\n#Programming\n#3D Plotting\n#Image Processing\n#Sounds\n#Videos\n#GUIs\n#Vectorization
MATLAB was designed to handle matrices, but one can do simple arithmetical operations:\n\n{{{>>5*3+6^2}}}\n{{{>>sin(pi/2)}}}\n\nYou can define a row vector by typing\n\n{{{>>[1 2 3]}}}\n\nYou can even name it \n\n{{{>>a=[1 2 3]}}}\n\nor, name it but not display the result\n\n{{{>>a=[1 2 3];}}}\n\nA row vector is obtained by inserting semicolons:\n\n{{{>>b = [1; 2; 3]}}}\n\nNow you can do the same operations on each element. Typing \n\n{{{>>2*a}}} \n\nwill produce [2 4 6]. But trying to square each element using \n\n{{{>>a^2}}} \n\nresults in an error. You can avoid this error by using the operator .^:\n\n{{{>>a.^2}}}\n\nSince {{{a}}} is a 1x3 row vector and {{{b}}} is a 3x1 column vector, we can compute \n\n{{{>>b*a}}} \n\nto obtain the outer product resulting in a 3x3 matrix or we can compute \n\n{{{>>a*b}}}\n\nto obtain the inner product. \n
[[Getting Started]]\n[[Arithmetic Operations]]\n[[Built-in Functions]]\n[[Matrices]]\n[[2D Plotting]]\n[[Scripts]]\n[[Functions]]
There are many 'Built-in Functions' in MATLAB. Some of the more commonly used functions are listed below. One can compute the value of a function appled to all elements of an array. For example, computing the exponential of each element of {{{A}}} is found by \n\n{{{>>exp(A)}}}\n\n\n|!Function |!Description |!Function|!Description |!Function |!Description |\n|acos |Inverse cosine |acosh |Inverse hyperbolic cosine|acot |Inverse cotangent|\n|acoth |Inverse hyperbolic cotangent |acsc |Inverse cosecant |acsh |Inverse hyperbolic cosecant|\n|asec |Inverse secant |asech |Inverse hyperbolic secant |asin |Inverse sine |\n|asinh |Inverse hyperbolic sine |atan |Inverse tangent |atanh |Inverse hyperbolic tangent|\n|atan2 |Four-quadrant inverse tangent|cos |Cosine|cosh |Hyperbolic cosine |\n|cot |Cotangent |coth |Hyperbolic cotangent |csc |Cosecant |\n|csch |Hyperbolic cosecant |hypot |Square root of sum of squares |sec |Secant |\n|sech |Hyperbolic secant |sin |Sine |sinh |Hyperbolic sine |\n|tan |Tangent |tanh |Hyperbolic tangent |exp |Exponential|\n|log |Natural logarithm |log2 |Base 2 logarithm|log10 |Common (base 10) logarithm|\n|sqrt |Square root|nthroot |Real nth root| abs |Absolute value|\n|angle |Phase angle |conj |Complex conjugate|i |Imaginary unit|\n|imag|Complex imaginary part|j|Imaginary unit| real|Complex real part|\n|sign |Signum|fix |Round towards zero |floor |Round towards minus infinity|\n|ceil |Round towards plus infinity |round |Round towards nearest integer |mod |Modulus after division|\n|rem |Remainder after division|factor |Prime factors |factorial |Factorial function |\n|gcd |Greatest common divisor |isprime |Determine if input is prime number|lcm |Least common multiple|\n|nchoosek |All combinations of N elements taken K at a time |perms |All possible permutations | primes|Generate list of prime numbers|\n|rat, rats |Rational fraction approximation| | | | |
Here are some things that did not work in the MATLAB Examples at this site:\n\n{{{plot(rvals,xresults,'b.','MarkerSize',1)}}} - Octave gives an error when setting Marker Size. Use {{{plot(rvals,xresults,'b.')}}} instead. First seen in ''bifurf''.
''The Logistic Map''\nThe typical example for 1D discrete maps is the logistic map, using the function //f(x)=rx(1-x)// defined for //x// in [0,1] and //r// in (0,4]. Here we show howto use MATLAB to study this map.\n\nOne can quickly display the orbits in the command window using the following lines: \n\n{{{>>f = @(x,r) r*x*(1-x)}}}\n{{{>>x(1)=0.1; for i=2:200, x(i)=f(x(i-1),3.5); end;}}}\n{{{>>plot(100:200,x(100:200),'o'), xlabel('n'), ylabel('x_n')}}}\n\nThis produces\n\n[img[Orbit of Logistic Map|logisticorbit.jpg]]\n\nA function can be written which displays a little more:\n\n{{{\nfunction y=logorbits(r,N)\n% This function computes an orbit of the logistic map.\n% r = parameter and N = number of elements of the orbit\n% It then displays the last 10 elements of the orbit\n% and a plot of the last 100 elements.\nf = @(x,r) r*x*(1-x);\nx(1)=0.1;\nfor i=2:N \n x(i)=f(x(i-1),r); \nend;\nplot(N-99:N,x(N-99:N),'o')\nxlabel('n')\nylabel('x_n')\ny=x(N-9:N);\n}}}\n\nIn the above figure we can see a period 4 orbit. We can plot //f(x)// and //f^^4^^(x))// with //y=x// to see that this is the case. \n\n{{{>>f=@(x,r) r*x.*(1-x)}}}\n{{{>>x=linspace(0,1,200);}}}\n{{{>>plot(x,x,x,f(x,3.5),x,f(f(f(f(x,3.5),3.5),3.5),3.5))}}}\n\n[img[Plot of f(x) and f^4(x)|period4.jpg]]\n\nComputation of iterations of //f(x)// can be complicated to type. So, we could define a function to handle such iterations.\n\n{{{\nfunction y=logisticN(x,r,N)\ny=x;\nfor i=1:N\n y=r*y.*(1-y);\nend\n}}}\n\nWhen we are interested in the stability of a given periodic, or fixed, point, we can just iterate starting at or near the point. If it is stable, the orbit will tend towards that point. Looking at how the orbits behave asymptotically as //r// is varied leads to what is call a bifurcation diagram. In the following code the orbits are followed for a long time and one needs only look at the last points. Graphing these vs //r// displays the patterns given in the figures. The second figure is a zoom on a section of the first diagram indicating the complexity at small scales.\n\n{{{\nfunction bifurf(ra,rb,xa,xb,maxiter)\n% Function evaluation is given by\n% bifurf(1,4,0,1,200)\n% bifurf(3.57,3.58,0.33,0.4,500000)\n%\n% modification of code in Matlab Guide, Higham and Higham, SIAM2005 \nr=linspace(ra,rb,1001)';\nic=0.2:0.3:1;\npts=50;\nrvals=repmat(r,length(ic),1);\nxresults=zeros(length(rvals),pts);\nx=kron(ic',ones(size(r)));\n\nfor k=2:maxiter\n x=rvals.*x.*(1-x);\nend\nfor k=1:pts\n x=rvals.*x.*(1-x);\n xresults(:,k)=x;\nend \n\nplot(rvals,xresults,'b.','MarkerSize',1)\naxis([ra,rb,xa,xb]);\nxlabel('r');\nylabel('x');\ntitle('Bifurcation for Logistic Map')\n}}}\n\n[img[Birfurcation Diagram for the Logistic Map|bifur.jpg]]\n\n[img[Birfurcation Diagram for the Logistic Map|bifur2.jpg]]\n\nAnother thing you might want to do is to find the periodic orbits. Recall that this means solving the fixed point equation //f^^4^^(x)=x.// To do this we define //g(x) = f^^4^^(x)-x.// Then we seek the zeros of //g(x)//. \n\nTo do this we first define //g// by typing \n\n{{{ >>g = @(x) logisticN(x,3.5,4)-x;}}}\n\nThis assumes that you have saved the function {{{logisticN}}}. To find the roots of //g(x)=0// you can use the function {{{fzero}}}. To see an example, type \n\n{{{>>help fzero}}}\n\nNot that you need to provide an interval in which you know there is a zero, or at least put in an approximation to the zero. So, type \n\n{{{>>fzero(g,0.4)}}}\n\nor\n\n{{{>>fzero(g,[0.3 0.4])}}}\n\nThe output is\n\n{{{ans = \n 0.3828}}}\n\nMATLAB displays only four digits by default, though it uses more digits that that in the background. To display more digits enter\n\n{{{>format long}}\n\nRecomputing the answer, you get\n\n{{{ans = \n 0.38281968301732}}}\n\nTo obtain the complete orbit, insert this result as //x~~0~~// into the logistic function several times. you can do this by typing\n\n{{{>>logisticN(0.38281968301732,3.5,1)}}}\n\nor, if it is the line right after computing the periodic point, you can type\n\n{{{>>logisticN(ans,3.5,1)}}}\n\nThe advantage of the last method is that you can use the up arrow and do this several times to generate the entire orbit. This gives the period four orbit {0.38281968301732, 0.82694070659144, 0.50088421030722, 0.87499726360246}.\n\nThere are other solutions. Recall that the fixed points of //f(x)// are also solutions. So, we know that //x = 0// and //x = 1 - 1/r = 1 - 1/3.5 = 0.7142857...// are solutions. Since 4 is divisible by 2, we expect a period two orbit. So, searching for intersections of //y = f^^4^^(x)// with //y = x//, we see that we have not captured the intersection for //x// in [0.4, 0.5]. Let's look for this point:\n\n{{{>> fzero(g,[0.4 0.5])}}}\n\nWe find //x = 0.42857142857143//. Following the above procedure for generating the orbit, we obtain \n{0.42857142857143, ,0.42857142857143, }. This is obviously the period two orbit we expected.\n \nWe can do all of this by placing our commands into a script for use at other times without doing all of the typing. Copy the folowing code and save as {{{periodic.m}}}. Then you need only only change //r// and //N// and the interval. Save the file and type {{{periodic}}} in the command window.\n{{{\n% Find periodic orbits for the the logitic map.\n%\n% N = period\n% r = parameter\n% [x0 x1] = search interval\n% y = set of points in orbit\n\n% Enter N and r\nN=4;\nr=3.5;\n\n% Plot y=f^N(x) and y=x.\nx=linspace(0,1,200);\nplot(x,x,x,logisticN(x,r,N))\n\n% Define f^N(x)-x = g =0.\ng = @(x) logisticN(x,r,N)-x;\n\n% Find root of g near x=0.4.\nX(1)= fzero(g,0.4);\n\n% Iterate logistic map to get orbit\nfor k=2:N\n X(k)=logisticN(X(k-1),r,1);\nend\n\n% Output the orbit\nformat long;\nX\nformat\n}}}\n\nNote that comments were added using the percent sign. You need only change the parameters of the search. \n\nThe next thing you might want to do is to apply these methods for other maps. In this case you need to replace {{{logisticN}}} with a new function. There are ways in MATLAB to generalize this for any map in order to simplify changes from one problem to the next. However, this would fall into more advanced MATLAB programming.
''Maps in the Plane''\n\nOne can now think about 2D maps, like the henon map. This is not the most efficient use of Matlab, but shows the basic aspects of programming.\n\nCreate a file called [[henon.m|../chaos/henon.m]]\n\n{{{\nN=500;\na=0.5;\nb=.3;\naxis([-2.5 2.5 -2.5 2.5]);\n[x0 y0 bb]=ginput(1); (Line for mouse input)\nhold\nwhile bb<2, (Left click to run, right click to end)\n plot(x0,y0);\n for i=1:N,\n x1=a-x0^2+b*y0;\n y1=x0;\n x0=x1;\n y0=y1;\n plot(x0,y0);\n [x0,y0]\n end\n [x0 y0,bb]=ginput(1);\nend\nhold;\n}}}\n\nAnother version of henon can do the iteration of the map for a set of initial conditions and plot the eigenvectors on the plot as well.\n\nCreate a file called [[henon2.m|../chaos/henon2.m]]\n\n{{{\na=1.4;\nb=.3;\nr=0.1;\nclf;\naxis([-2.5 2.5 -2 2]);\ntitle(['Henon Map ' num2str(N) ' Iterations']);\nhold on\nfor i=0:200\n x0 = .883896+r*cos(i*2*pi/200);\n y0 = .883896+r*sin(i*2*pi/200);\n plot(x0,y0);\n for i=1:N,\n x1=a-x0^2+b*y0;\n y1=x0;\n x0=x1;\n y0=y1;\n if i==N\n plot(x0,y0);\n end\n end\nend\nline([.883896+0,.883896+.155946],[.883896+0,.883896+1]);\nline([.883896+0,.883896-1.9237],[.883896+0,.883896+1]);\nhold off;\n}}}\n\n!!Template for 2D maps:\n\nWe can do a general 2D map like the 1D map above. Here is a map called the standard map. Put it into a function file called [[map2d.m|../chaos/map2d.m]]\n\n{{{\nfunction [x1,y1]=map2d(x0,y0)\nk=0.7;\ny1=y0-k/(2*pi)*sin(2*pi*x0);\nx1=y1+x0;\ny1=mod(y1,1);\nx1=mod(x1,1);\n}}}\n\nWe then write a routine to iterate this map and plot the points as it iterates. Save this under a name like [[iter2d.m|../chaos/iter2d.m]] and type iter2d in the Command Window.\n\n{{{\nx(1)=0.0;\ny(1)=0.3;\naxis([0 1 0 1]);\nhold on\np = plot(x(1),y(1),'.','EraseMode','none','MarkerSize',2); \nfor n=2:1000\n [x(n),y(n)]=map2d(x(n-1),y(n-1));\n set(p,'XData',x(n),'YData',y(n))\n drawnow\nend\nhold off\n}}}\n \nA modification of this routine will allow us to plot orbits for various initial conditions. We add ginput and a while loop.\n\n{{{\naxis([0 1 0 1]);\n[x(1) y(1) bb]=ginput(1);\nhold on\np = plot(x(1),y(1),'.','EraseMode','none','MarkerSize',2); \nwhile bb<2,\n for n=2:500\n [x(n),y(n)]=map2d(x(n-1),y(n-1));\n set(p,'XData',x(n),'YData',y(n))\n drawnow\n end\n [x(1) y(1) bb]=ginput(1);\nend\nhold off\n}}}\n\nAlso see [[gingermap.m|../chaos/gingermap.m]], [[ginger.m|../chaos/ginger.m]] for the gingerbread man map.
MATLAB comes with common [[Built-in Functions]] and many other specialized functions in a variety of areas. Some of these can be found by looking at MATLAB's ''Help'' section. If these do not satisfy your needs, you can design your own functions. Functions take some input, in the form of scalars, vectors, or matrices, and can output specific results. \n\nAs a simple example, lets say that we want the area of a triangle give the three sides, //a//, //b//, //c//. As you know, this can be solved using [[Heron's Formula|http://mathworld.wolfram.com/HeronsFormula.html]]: //A = (s(s-a)(s-b)(s-c))^^1/2^^ //, where // s = (a+b+c)/2//.\n\nThe input variables are the sides and the output is the area. Open up the editor and type in the following lines and save as heron.m:\n\n{{{\nfunction f=heron(a,b,c)\n% This function computes the area of a triangle with sides a,b,c\n\ns=(a+b+c)/2;\nf=sqrt(s*(s-a)*(s-b)*(s-c));\n}}}\n\nIn the command window you can call this function:\n\n{{{>>heron(3,4,5)}}}\n\nIt produces an answer of {{{6}}}. Since this is a known triangle, it is easy to see that the function works. Now you can try the function with several other triangles. You can also call this function from other scripts and functions provided heron.m is in the path. \n\nIn the above function we have added a comment statement. This comment will appear if you now type in \n\n{{{>>help heron}}}\n\nAlso, we were careful to add semicolons after each line so that we do not get extraneous output displayed in the command window.\n\nIf we type in any of the variables {{{a}}}, {{{b}}}, {{{c}}}, {{{s}}}, we will not get any output. These are what are called local variables. One can get a bit more sophisticated by declaring global variables known to several functions and scripts. One can also pass parameters, input matrices, or output matrices. But the simplest form of functions has been provided by this example. \n\nSometimes one wants to define a one line function and creating a file might seem like a lot to do. In more recent version of MATLAB the //inline// object has been introduced. One can define a function, like //f(x)=exp(-x^^2^^) //, simply by typing in the command window or in some script:\n\n{{{>>f = inline(exp(-x^2);}}} \n\nNow evaluate //f(0.1)//. However, one cannot use this function on a vector //x// since the //x^^2^^// would return an error. In such cases you need to add a dot:\n\n{{{>>f = inline(exp(-x.^2);}}}\n\nThen you could do the following:\n\n{{{>>x=[0.1 0.2 0.3];}}}\n{{{>>f(x)}}} \n\nto obtain {{{[ 0.9900 0.9608 0.9139]}}}.\n\nOne final way to handle one line functions is by using //anonymous functions//, which are a recent addition to MATLAB. This can be done for the above example:\n\n{{{>>f = @(x) exp(-x^2);}}} \n{{{>>f(0.1)}}}\n\nOne can define functions of several variables with a parameter previously defined in the script or work session:\n\n{{{>>a = 6.0}}}\n{{{>>f = @(x,y) x^2 + a*y^2}}}\n{{{>>f(2,3)}}}
A GUI is a Graphical User Interface. It is an advanced feature which allows users to design a control panel to run their applications with buttons, sliders, etc to manipulate parameters in their routines.
We will be using MATLAB 7 for MS Windows. When you launch the program you will get a main window with several subwindows. The subwindow on the right is the Command Window. You can type simple commands in this window after >> In this window you can try simple commands. Note that MATLAB is case sensitive. \n\nOne way to learn about MATLAB is to try the demos. Type demo and hit RETURN. Then under the MATLAB folder on the left you can select different demos. \n\nIf you want HELP, you can access topics under the menu item or type help topic in the command line. \nYou can recall previously types lines using the up and down arrows. You can quit by typing exit or quit, or just exit from the menu.
Iterated Function Systems\n\nIterated functions systems are just another example of iteration. As an example, here is the IFS code for generating a fern and watching the evolution of the plot. Again, create a file with a name like fern.m.\n\n{{{\n y=[0,0];\n p = plot(y(1),y(2),'.','EraseMode','none','MarkerSize',2); \n N=5000;\n xa=-5;\n xb=5;\n ya=0;\n yb=12;\n axis([xa xb ya yb]);\n hold on\n a=[0,.85,0.2,-0.15];\n b=[0,.04,-.26,.28];\n c=[0,-.04,.23,.26];\n d=[.16,.85,.22,.24];\n e=[0,0,0,0];\n f=[0,1.6,1.6,.44];\n Pr=[.01,.85,.07,.07];\n x0 = 0;\n y0 = 0;\n for i=1:N\n r = rand;\n if r < Pr(1) \n x1 = a(1) * x0 + b(1) * y0 + e(1);\n y1 = c(1) * x0 + d(1) * y0 + f(1);\n elseif r < Pr(1) + Pr(2) \n x1 = a(2) * x0 + b(2) * y0 + e(2);\n y1 = c(2) * x0 + d(2) * y0 + f(2);\n elseif r < Pr(1) + Pr(2)+Pr(3) \n x1 = a(3) * x0 + b(3) * y0 + e(3);\n y1 = c(3) * x0 + d(3) * y0 + f(3);\n else\n x1 = a(4) * x0 + b(4) * y0 + e(4);\n y1 = c(4) * x0 + d(4) * y0 + f(4);\n end \n x0 = x1;\n y0 = y1;\n set(p,'XData',x0,'YData',y0)\n drawnow\n end\n hold off\n}}}\n\nA slightly more compact form of this code is given below. Notice the use of matrices and commented lines.\n\n{{{\n x=[0;0];\n p = plot(x(1),x(2),'.','EraseMode','none','MarkerSize',2); \n N=5000;\n axis([-5 5 0 12]);\n hold on\n % Enter each affine transformation Ax+T\n A1=[0 0; 0 .16];\n T1=[0;0];\n A2=[.85 .04; -.04 .85];\n T2=[0;1.6];\n A3=[0.2 -.26; .23 .22];\n T3=[0;1.6];\n A4=[-.15 .28; .26 .24];\n T4=[0;.44];\n % Enter the probabilities\n Pr=[.01,.85,.07,.07];\n % The iteration\n for i=1:N\n r = rand;\n if r < Pr(1) \n y = A1*x+T1; \n elseif r < Pr(1) + Pr(2) \n y = A2*x+T2; \n elseif r < Pr(1) + Pr(2)+Pr(3) \n y = A3*x+T3; \n else\n y = A4*x+T4; \n end \n x=y;\n set(p,'XData',y(1),'YData',y(2))\n drawnow\n end\n hold off\n}}}
''MATLAB'' is a powerful envronment for doiing numerical computations. However, its one major drawback is that it is expensive. A student can obtain the [[Student Version on MATLAB|http://www.mathworks.com/academia/student_version/]]. Other options are to get one of the [[MATLAB clones|http://www.dspguru.com/sw/opendsp/mathclo2.htm]]: [[Octave|http://www.gnu.org/software/octave/]], [[SciLab|http://www.scilab.org/]], and [[RLab|http://rlab.sourceforge.net/]]. \n\nHere are some notes on\n[[Octave]] \n[[SciLab]]
[[What is MATLAB?]]\n[[MATLAB Clones]]\n-------------------------------\n[[Basic MATLAB]]\n[[Advanced MATLAB]]\n------------------------------\n''Topics''\n[[Linear Algebra]]\n[[Discrete Dynamics I]]\n[[Discrete Dynamics II]]\n[[IFS]]\n[[Numerical Analysis]]\n[[ODE Examples]]\n[[PDE Examples]]\n[[Signal Analysis]]\n-------------------------------\n[[Other Tutorials]]\n[[About]]\n-----------------------\n{{small{\n<<search>>\n<<closeAll>>\n<<permaview>>\n<<newTiddler>>\n<<saveChanges>>}}}\n---------------------------\n{{smaller{[[Site created and maintained by\n Dr. R. Herman.|http://people.uncw.edu/herman]]}}}
MATLAB's power is that one can work with matrices. In the section [[Arithmetic Operations]] we have described how to enter row and column vectors, which are 1x//n// and //n//x1 matrices, respectively. In order to build bigger matrices, one can type in the matrices row by row. So,\n\n{{{>>[1 2 3; 4 5 6; 7 8 9]}}} \n\nwill produce the matrix \n\n{{{ans=}}} \n\n{{{ 1 2 3}}}\n{{{ 4 5 6}}}\n{{{ 7 8 9}}}\n\nWriting \n\n{{{>>A=[1 2 3; 4 5 6; 7 8 9]}}} \n\nwe can call up a particular element using {{{A(2,3)}}} to extract the second row and third column, {{{6}}}. you can also change and element by entering \n\n{{{>>A(2,3)=10;}}}\n{{{>>A}}}\n\nto obtain\n\n{{{A=}}} \n\n{{{ 1 2 3}}}\n{{{ 4 5 10}}}\n{{{ 7 8 9}}}\n\nYou can also extract rows and columns. \n\n{{{>>B=A(1,:)}}}\n{{{>>C=A(:,1}}}\n\n{{{B}}} is the first row and {{{C}}} is the first column. \n\nTyping\n\n{{{>>A(1,:)=2*B}}} yields \n\n{{{ans=}}} \n\n{{{ 2 4 6}}}\n{{{ 4 5 10}}}\n{{{ 7 8 9}}}\n\nThere are several special matrices: {{{eye(N)}}} is the //N//x//N// identity matrix, {{{zeros(N,M)}}} produces and //N//x//M// matrix of zeros, {{{ones(N,M)}}} produces and //N//x//M// matrix of ones, and {{{diag(a)}}} produces a diagonal matrix with vector {{{a}}} along the main diagonal.
!!Systems of ODEs\n\n We will see shortly that systems of ODEs describe dynasmical systems and have all of the wonderful types of orbits we have been studying. You might want to numerically solve such systems. Here is an example of solving a system of two first order differential equations. The system is\n\n x'(t) = y(t)\n y'(t) = - x(t).\n\nFirst create a function file, like odefunc.m as shown below.\n\n{{{\n function dy = odefunc(t,y)\n dy = zeros(2,1); % a column vector\n dy(1) = y(2);\n dy(2) = -y(1);\n}}}\n\nNow you can use one of many ode solvers provided by Matlab, or write your own. We will use ode45, Type on the Command line\n\n [T,Y] = ode45('odefunc',[0 1000],[2 0]);\n\n To get a phase portrait, type plot(Y(:,1),Y(:,2)), or to get a solution vs time, type plot(T,Y(:,1)).\n\n To plot a family of solutions of a first order ODE, one can use the following set of files. The first file, which can be saved as ode1d.m, calls the differential equation from the function file odefun1.m. Clicking on the graph produces a solution curve for the solution passing through the point (t0,y(t0).\n\n{{{\n clf\n axis([0 100 -1 2]);\n [t0 x0 bb]=ginput(1); \n hold on\n while bb<2, \n [T,Y] = ode45('odefun1d',[t0 100],x0); \n plot(T,Y), \n [t0 x0,bb]=ginput(1);\n end\n hold off;\n\n function dy = odefun1d(t,y)\n dy = 0; \n dy = y*(1-y);\n}}}\n\nHere is a routine for generating a variety of initial condtions for a system of differential equations specified in a file odefunc.m. Both the phase portrait and the soltion are plotted in two windows. In this example a system needs to be input. An example of a linear systems follows.\n\n{{{\n % Routine to plot phase portrait and solutions to system odefunc\n % Plot window h1 is the solution window\n % Plot window h2 is the portrait window\n % Note the new functions for plotting windows:\n % subplot, title\n % \n h1=subplot('Position',[0.1 0.03 0.8 0.18]);\n axis([0 100 -1 1]); % axis for solution\n title(['Solution vs Time']);\n subplot(h1)\n hold on\n h2=subplot('Position',[0.15 0.27 0.5 0.7]);\n axis([-5 5 -5 5]); % axis for portrait\n title(['y(2) vs y(1)']);\n subplot(h2)\n [x0 y0 bb]=ginput(1); \n hold on\n while bb<2, \n [T,Y] = ode45('odefunc',[0 100],[x0 y0]); \n subplot(h2)\n plot(Y(:,1),Y(:,2)), % plot portrait\n subplot(h1)\n plot(T,Y(:,1)), % plot solution\n subplot(h2)\n [x0 y0,bb]=ginput(1);\n end\n hold off;\n subplot(h1)\n hold off\n}}}\n\nHere is the set of function for a linear system of ODEs.\n\n{{{\n function dy = odefunc(t,y)\n dy = zeros(2,1); % a column vector\n a=-1.0;\n b=-1.0;\n c=1.0;\n d=0.5;\n dy(1) = a*y(1)+b*y(2);\n dy(2) = c*y(1)+d*y(2);\n}}}