October 2009 – Numerical Methods Guy

How do I numerically solve an ODE in MATLAB?

The other day a student came to ask me for help in solving a second order ordinary differential equation using the ode45 routine of MATLAB. To use ode45, one needs to be familiar with how the inputs are required by MATLAB. The understanding of these inputs is important to use ode45 successfully in problems that are more complex than solving a second order ODE.

The ordinary differential equation was
2y”+3y’+5y=7 exp(-x), y(0)=11, dy/dx(0)=13
This has to put in the state variable form by reducing it by using
y’=z
That gives
y’=z with the corresponding initial conditions as y(0)=11
Then
2y”+3y’+5y=7 exp(-x)
reduces to
2z’ + 3z+5y=7exp(-x)
z’ =(7exp(-x)-3z-5y)/2 with the corresponding initial conditions as z(0)=13

So as needed by MATLAB, call y as y(1) and z as y(2)
dy(1)=y(2), y(1) at x=0 is 11
dy(2)=(7exp(-x)-3y(2)-5y(1))/2, y(2) at x=0 is 13

These equations are now put in a MATLAB function we call odestate.m
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=(7*exp(-x)-3*y(2)-5*y(1))/2;

To solve the ODE, the
The inputs are
1) the function odestate
2) The outputs are required between x=0 and x=17,
hence entered as [0 17]
3) The initial conditions are y(0)=11 and dy/dx(0)=13,
hence entered as [11 13]

The outputs are
1) X= array of x values between 0 and 17
2) Y= matrix of 2 columns;
first column is the y(x)
second column is dy/dx(x)
The MATLAB code then is
[X,Y]=ode45(@odestate,[0 17],[11 13]);

Click the links for the MATLAB mfiles for the function odestate.m and the ODE solver odetest.m

_________________________________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://nm.mathforcollege.com/videos and http://www.youtube.com/numericalmethodsguy

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you

The continue statement in MATLAB

The continue statement in MATLAB is used to pass control to the next iteration in for and while statements. Let’s suppose someone wants to find and print the value of k^2-50 for all integers in [-10,10] domain. The mfile for that is given below.

% For integers k=-10,-9,….,9,10,
% the function k^2-50 will take positive as
% well as negative values.
%For example, for k=-9, k^2-50=31; for k=1,
% k^2-50=-49; for k=8, k^2-50=14.
% The loop below will calculate values of k^2-50 for
% all values of requested k.
for k=-10:1:10
val=k^2-50;
fprintf(‘\n k=%g val=%g’,k,val)
end

___________________________________________________________

Let’s suppose now you are asked to calculate and print value of k^2-50 for all integers in [-10,10] domain but only if (k^2-50) is positive.

% The loop below will calculate and print values of k^2-50
% for all values of the requested k
% for which k^2-50 is positive.
for k=-10:1:10
    if (k^2-50<0)
        continue;
    end
    val=k^2-50;
    fprintf(‘\n k=%g val=%g’,k,val)
end

____________________________________________________________

Can you do what you did above using the while statement. Yes, the MATLAB code is given below.

% Equivalent in while
% The loop below will calculate values of k^2-50
% for all values of the requested k for which k^2-50 is positive.
k=-10;
while (k<=10)
    if (k^2-50>0)
        val=k^2-50;
        fprintf(‘\n k=%g val=%g’,k,val)
    end
k=k+1;
end

_______________________________________________________

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

The break statement in MATLAB

The break statement in MATLAB is used to break out of a loop – a for or while statement, that is, it terminates the execution of the loop. Let’s suppose someone wants to find the value of k^2-50 for all integers in [-10,10] domain. The mfile for that is given below.

% For integers k=-10,-9,….,9,10,
% the function k^2-50 will take positive as
% well as negative values.
%For example, for k=-9, k^2-50=31; for k=1,
% k^2-50=-49; for k=8, k^2-50=14.
% The loop below will calculate values of k^2-50 for all values of requested k.
for k=-10:1:10
val=k^2-50;
end

___________________________________________________________

Let’s suppose now you are asked to calculate value of k^2-50 for all integers in [-10,10] domain but only until k^2-50 becomes negative.

% The loop below will calculate values of k^2-50
% for all values of the requested k
% until it turns negative
for k=-10:1:10
    if (k^2-50<0)
        break;
    end
    val=k^2-50;
    fprintf(‘\n k=%g val=%g’,k,val)
end

____________________________________________________________

Can you do what you did above using the while statement. Yes, the MATLAB code is given below.

% Equivalent in while
% The loop below will calculate values of k^2-50
% for all values of the requested k until it turns negative
k=-10;
while (k<=10) & (k^2-50>0)
    val=k^2-50;
    fprintf(‘\n k=%g val=%g’,k,val)
    k=k+1;
end

_______________________________________________________

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.