Introduction to programming in Matlab 2

by Kristian Sandberg
Department of Applied Mathematics
University of Colorado

In this worksheet we will continue to learn Matlab programming. The main goal is that you by completing the worksheet can write your own Runge Kutta 4 ODE-solver. The worksheet centers around a few examples, and it is important to implement these examples, run them and carefully compare the output with the code.

Functions in Matlab

You have already encountered functions that are built in to Matlab. The sin() is a function that takes an argument and returns an output. As you start to write your own programs in Matlab you will need to create your own functions that takes one or more arguments and does something with these arguments. For example, when you write an ODE-solver you will probably have a loop calling "the right-hand side" of your ODE (i.e., f(t,x) in x'=f(t,x)) thousands of time. You could type this function explicitly in your loop, but you probably want your ODE-solver to be able to handle more than one ODE and then it will convenient to write your right-hand side of the ODE as a function. Another example when you want to use functions in Matlab is if you want to make a graphing program. This will be illustrated in the example below, but let's first look at the syntax for a function in Matlab.

function y=function name(argument list)

commands

The code above must be written in a separate m-file! The name of the file should coincide with the name of the function. Remember to save the file with a .m after the file name. The syntax above is illustrated by an example. The example defines a function in a function file named f.m and is saved in the same directory as the file functionplotter.m.

Example 1.

% Filename: f.m function y = f(t) y = exp(-t)*sin(3*t);

% Filename: functionplotter.m clear for k=1:1000 t(k)=0.01*k; y(k)=f(t(k)); end plot(t,y)

You must type and save the file f.m before you can use the code in functionplotter.m. Once you have typed and saved both these files and saved them in the same directory, you run the program in the Matlab command window by typing functionplotter. This code will call the function file f.m each time it needs to compute the function you specified in the file f.m. Note that it is not necessary to have a clear command in the function file f.m. Make sure you understand how this example works. It will help you to write the ODE-solver later.

Exercise 1.
Make a plot of the function exp(t)*cos(5*t) for t between 0 and 10. (Hint: Do you need to change the file functionplotter.m at all?)

It is possible for functions to have more than one argument which the following "stupid" example illustrates.

Example 2.

% Filename: g.m (remember, file name must match % with the name of the function) function y = g(t,x1,x2) y = sin(t+x1+x2);

% Filename: test.m clear function_value = g(1,2,3)

The example above shows how you can use functions of three variables.

How to write large programs

When one has to write a "large" program the first time, it is easy to feel overwhelmed by the task. You ask yourself, "where should I start?". No matter how experienced programmer you are, you need to find a way to organize your work. Here are a few advice how you can organize your work.

Before you begin coding ask yourself (and write down your answer on a piece of paper) the following questions:

What is my input? What should I name these variables? What types are my inputs? (That is, are my inputs files, variables, vectors...?)
What is my output? (For example, is it a number, a plot or something else.)
What are my main procedures in my program? That is, what are the main steps my program needs to perform? (Don't answer in details and commands at this stage. Try to specify your answers to these questions in words.)

One way of answering these questions, is in a form of diagram. The diagram for an Euler forward ODE solver may look something like shown in Diagram 1.

Once the above questions are answered, we observe that the most tricky part might be the problem of iterating Eulers method to find a solution. Therefore we should try to specify this step. One way of specifying what we need to do in this step is to extend our first diagram by adding information of how to do the loop. This can be experessed symbolically as in Diagram 2.

Now we have specified what our program has to do so now we are in position to start coding. In this case, I would recommend you to start with the function file. Once this one is written, start typing your main code, but make sure to check that each (small) step works, before you proceed to the next step of the code! Here is an example of the code for an ODE solver. In case you feel comfortable with programming you can already now start doing Exercise 2 below, rather than typing in the code of Example 2. The following example solves the equation x'(t)=-x which you can solve analytically. The reason for choosing such an easy function is that you know what the solution should look like and therfore you can check that your solver works as expected before you try more tricky equations. Below I will use the % frequently to add comments to the code. I will also use the command ceil( ) which take an argument and rounds it to the nearest greater integer. Please compare the code to Diagram 2.

Example 3.

% Filename: f.m (name must match with the name of the function) function y = f(x) y = -x;

% Filename: ef.m clear t0 = 0; % Start time tfinal = 10; % Final time x0 = 1; % Initial value (x(t0)) h = 0.05; % Time step N = ceil( (tfinal-t0)/h ); % Computes the number of loop steps t(1) = t0; % t( ) will be a vector with the time grid x(1) = x0; % x( ) will be a vector with the "solution" for each time value for k=1:N t(k+1) = h*k; % Fills the time grid vector x(k+1) = x(k)+h*f(x(k)); % Euler forward iteration end plot(t,x)

Exercise2.
Write your own RK4 solver for autonomous first order ODEs. Test your solver on the following IVP: x'=x, x(0)=1. If it works, try it on some more tricky equation. The RK4 iteration scheme for autonomous first order ODEs is
k1=f(x(k))
k2=f(x(k)+h*k1/2)
k3=f(x(k)+h*k2/2)
k4=f(x(k)+h*k3)
x(k+1)=x(k)+h*(k1+2*k2+2*k3+k4)/6
If you find this exercise too difficult, you can work on the examples and exercises given in this and the other worksheets and maybe try to extend/improve them. Do not feel pressure to finish the RK4-solver but if you like programming this is a good exercise.