How do I solve an initial value ODE in MATLAB?

Many students ask me how do I do this or that in MATLAB.  So I thought why not have a small series of my next few blogs do that.  In this blog, I show you how to solve an initial value ordinary differential equation.

• The MATLAB program link is here.
• The HTML version of the MATLAB program is here.
• DO NOT COPY AND PASTE THE PROGRAM BELOW BECAUSE THE SINGLE QUOTES DO NOT TRANSLATE TO THE CORRECT SINGLE QUOTES IN MATLAB EDITOR.  DOWNLOAD THE MATLAB PROGRAM INSTEAD

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I solve an initial value ordinary differential equation?

%% SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/ode_initial.m;
% Last Revised : May 14, 2009;
% Abstract: This program shows you how to solve an
%           initial value ordinary differential equation.
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to solve’)
disp(‘   an initial value ordinary differential equation’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/ode_initial.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   May 14, 2009’)
disp(‘ ‘)

%% INPUTS
% Solve the ordinary differential equation 3y”+5y’+7y=11exp(-x)
% Define x as a symbol
syms x
%The ODE
ode_eqn=’3*D2y+5*Dy+7*y=11*exp(-13*x)’;
% The initial conditions
iv_1=’Dy(0)=17′;
iv_2=’y(0)=19′;
% The value at which y is sought at
xval=23.0;
%% DISPLAYING INPUTS

disp(‘INPUTS’)
func=[‘  The ODE to be solved is ‘ ode_eqn];
disp(func)
iv_explain=[‘  The initial conditions are ‘ iv_1 ‘    ‘ iv_2];
disp(iv_explain)
fprintf(‘  The value of y is sought at x=%g’,xval)
disp(‘  ‘)

%% THE CODE

% Finding the solution of the ordinary differential equation
soln=dsolve(ode_eqn,iv_1,iv_2,’x’);
% vpa below uses variable-precision arithmetic (VPA) to compute each
% element of soln to 5 decimal digits of accuracy
soln=vpa(soln,5);

%% DISPLAYING OUTPUTS
disp(‘  ‘)
disp(‘OUTPUTS’)
output=[‘  The solution to the ODE is ‘ char(soln)];
disp(output)
value=subs(soln,x,xval);
fprintf(‘  The value of y at x=%g is %g’,xval,value)
disp(‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos

How do I solve a nonlinear equation that needs to be setup in MATLAB?

Many students ask me how do I do this or that in MATLAB. So I thought why not have a small series of my next few blogs do that. In this blog, I show you how to solve a nonlinear equation that needs to be set up.

For example to find the depth ‘x’ to which a ball is floating in water is based on the following cubic equation
4*R^3*S=3*x^2*(R-x/3)
where
S= specific gravity of ball
So how do we set this up if S and R are input values?

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I solve a nonlinear equation if I need to set it up?

%% SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/integration.m;
% Last Revised : March 28, 2009;
% Abstract: This program shows you how to solve a nonlinear equation
% that needs to set up as opposed that is just given to you.
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to solve’)
disp(‘   a nonlinear equation that needs to be setup’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/nonlinearequation.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   April 17, 2009’)
disp(‘ ‘)

%% INPUTS
% Solve the nonlinear equation where you need to set up the equation
% For example to find the depth ‘x’ to which a ball is floating in water
% is based on the following cubic equation
% 4*R^3*S=3*x^2*(R-x/3)
% S= specific gravity of ball
% So how do we set this up if S and R are input values

S=0.6
R=0.055
%% DISPLAYING INPUTS
disp(‘INPUTS’)
func=[‘  The equation to be solved is 4*R^3*S=3*x^2*(R-x/3)’];
disp(func)
disp(‘  ‘)

%% THE CODE
% Define x as a symbol
syms x
% Setting up the equation
C1=4*R^3*S
C2=3
f=[num2str(C1) ‘-3*x^2*(‘ num2str(R) ‘-x/3)’]
% Finding the solution of the nonlinear equation
soln=solve(f,x);
solnvalue=double(soln);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
for i=1:1:length(solnvalue)
fprintf(‘\nThe solution# %g is %g’,i,solnvalue(i))
end
disp(‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos

How do I solve a nonlinear equation in MATLAB?

Many students ask me how do I do this or that in MATLAB.  So I thought why not have a small series of my next few blogs do that.  In this blog, I show you how to solve a nonlinear equation.

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I solve a nonlinear equation?

%% SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/integration.m;
% Last Revised : March 28, 2009;
% Abstract: This program shows you how to solve a nonlinear equation.
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to solve’)
disp(‘   a nonlinear equation’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/nonlinearequation.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   April 11, 2009’)
disp(‘ ‘)

%% INPUTS
% Solve the nonlinear equation x^3-15*x^2+47*x-33=0
% Define x as a symbol
syms x
% Assigning the fleft hand side o the equation f(x)=0
f=x^3-15*x^2+47*x-33;
%% DISPLAYING INPUTS

disp(‘INPUTS’)
func=[‘  The equation to be solved is ‘ char(f), ‘=0’];
disp(func)
disp(‘  ‘)

%% THE CODE

% Finding the solution of the nonlinear equation
soln=solve(f,x);
solnvalue=double(soln);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
for i=1:1:length(solnvalue)
fprintf(‘\nThe solution# %g is %g’,i,solnvalue(i))
end
disp(‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos

How do I integrate a discrete function in MATLAB?

Many students ask me how do I do this or that in MATLAB.  So I thought why not have a small series of my next few blogs do that.  In this blog, I show you how to integrate a discrete function.

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

_____________________________________________________

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I integrate a discrete function?  Three cases of data are
% discussed.

%% SUMMARY

% Language : MATLAB 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/integrationdiscrete.m;
% Last Revised : April 3, 2009;
% Abstract: This program shows you how to integrate a given discrete function.

clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to integrate’)
disp(‘   a discrete function’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/integrationdiscrete.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   April 3, 2009’)
disp(‘ ‘)

%% CASE 1

%% INPUTS

% Integrate the discrete function y from x=1 to 6.5
% with y vs x data given as (1,2), (2,7), (4,16), (6.5,18)
% Defining the x-array
x=[1  2  4  6.5];
% Defining the y-array
y=[2  7  16  18];

%% DISPLAYING INPUTS
disp(‘____________________________________’)
disp(‘CASE#1’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION MATCH x(1) AND x(LAST)’)
disp(‘ ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,x(1))
fprintf(‘\n  Upper limit of integration, b= %g’,x(length(x)))
disp(‘ ‘)

%% THE CODE

intvalue=trapz(x,y);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)
disp(‘___________________________________________’)

%% CASE 2

%% INPUTS

% Integrate the discrete function y from x=3 to 6
% with y vs x data given as (1,2), (2,7), (4,16), (6.5,18)
% Defining the x-array
x=[1  2  4  6.5];
% Defining the y-array
y=[2  7  16  18];
% Lower limit of integration, a
a=3;
% Upper limit of integration, b
b=6;
%% DISPLAYING INPUTS

disp(‘CASE#2’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST)’)
disp(‘  ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
% Choose how many divisions you want for splining from a to b
n=1000;
fprintf(‘\n  Number of subdivisions used for splining = %g’,n)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE

xx=a:(b-a)/n:b;
% Using spline to approximate the curve from x(1) to x(last)
yy=spline(x,y,xx);
intvalue=trapz(xx,yy);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)
disp(‘___________________________________________’)
%% CASE 3

%% INPUTS

% Integrate the discrete function y from x=1 to 6.5
% with y vs x data given as (1,2), (4,16), (2,7), (6.5,18)
% The x-data is not in ascending order
% Defining the x-array
x=[1  4   2 6.5];
% Defining the y-array
y=[2  16  7 18];
% Lower limit of integration, a
a=3;
% Upper limit of integration, b
b=6;
%% DISPLAYING INPUTS

disp(‘CASE#3’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST) ‘)
disp(‘AND X-DATA IS NOT IN ASCENDING OR DESCENDING ORDER’)
disp(‘   ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
% Choose how many divisions you want for splining from a to b
n=1000;
fprintf(‘\n  Number of subdivisions used for splining = %g’,n)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE
[x,so] = sort(x); % so is the sort order
y = y(so); % y data is now in same order as x data
xx=a:(b-a)/n:b;
% Using spline to approximate the curve from x(1) to x(last)
yy=spline(x,y,xx);
intvalue=trapz(xx,yy);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)

____________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos

How do I integrate a continuous function in MATLAB

Many students ask me how do I do this or that in MATLAB.  So I thought why not have a small series of my next few blogs do that.  In this blog I show you how to integrate a continuous function.

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

___________________________________________

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I integrate a continuous function?

%% SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/integration.m;
% Last Revised : March 28, 2009;
% Abstract: This program shows you how to integrate a given function.
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to integrate’)
disp(‘   a given function ‘)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/integration.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   March 29, 2009’)
disp(‘ ‘)

%% INPUTS

% Integrate exp(x)*sin(3*x) from x=2.0 to 8.7
% Define x as a symbol
syms x
% Assigning the function to be differentiated
y=exp(x)*sin(3*x);
% Assigning the lower limit
a=2.0;
% Assigning the upper limit
b=8.7;

%% DISPLAYING INPUTS

disp(‘INPUTS’)
func=[‘  The function is to be integrated is ‘ char(y)];
disp(func)
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE

% Finding the integral using the int command
% Argument 1 is the function to be integrated
% Argument 2 is the variable with respect to which the
%    function is to be integrated – the dummy variable
% Argument 3 is the lower limit of integration
% Argument 4 is the upper imit of integration
intvalue=int(y,x,a,b);
intvalue=double(intvalue);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)

_________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos

MATLAB code for the efficient automatic integrator

In the previous post, we discussed why doubling the number of segments in the automatic integrator based on multiple-segment trapezoidal rule is more efficient than increasing the number of segments one at a time. But this advantage involves having to store the individual function values from previous calculations and then having to retrieve them properly. This drawback was circumvented very efficiently by using the formula derived in another previous post where there is no need to store individual function values.

The matlab file for finding a definite integral by directly using the multiple segment trapezoidal rule from this post is given here (matlab file, html file), while the matlab file that uses the more efficient formula from this post is given here (matlab file, html file).  Here are the inputs to the programs.

% a = Lower limit of integration
% b = Upper limit of integration
%  nmax = Maximum number of segments
% tolerance = pre-specified tolerance in percentage
% f = inline function as integrand

a=5.3;
b=10.7;
nmax=200000;
tolerance=0.000005;
f=inline(‘exp(x)*sin(2*x)’)

We ran both the program on a PC and found that the more efficient algorithm (51 seconds) ran in half the time as the other one (82 seconds).  This is expected, as only n function evaluations are made for 2n-segments rule with the efficient formula, while 2n+1 functions evaluations are made for the original formula.

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

Is a square matrix diagonal or not?

A square matrix A is diagonal if all the elements on the off-diagonal are zero. That is, A(i,j)=0 for i~=j.

In this posting, I show a MATLAB program that finds whether a square matrix is diagonal by using three different methods. These are academic ways to reinforce programming skills in a student.

The MATLAB program can be downloaded as a Mfile (better to download it, as single quotes from the web-post do not translate correctly with the MATLAB editor). The html file showing the mfile and the command window output is also available.

%% IS A GIVEN SQUARE MATRIX A DIAGONAL MATRIX?
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : November 15, 2008
% Abstract: This program shows you three ways of finding out
% if a square matrix is a diagonal matrix. A square matrix is
% diagonal if all the off-diagonal elements are zero, that is
% A(i,j)=0 for i~=j.
clc
clear all
disp(‘This program shows you three ways of finding out’)
disp(‘if a square matrix is a diagonal matrix.’)
disp(‘A square matrix is diagonal if all the off-diagonal’)
disp(‘elements are zero, that is A(i,j)=0 for i~=j.’)
disp(‘ ‘)
%% INPUTS
% The square matrix
A=[1 0 0 0;0 3.4 0 0; 0 0 -4.5 0;0 0 0 0];
disp (‘INPUTS’)
disp(‘Here is the square matrix’)
A
disp(‘ ‘)

%% FIRST SOLUTION
% This is based on counting the number of zeros on
% off the diagonal. If this count is n^2-n then it
% is a diagonal matrix, otherwise it is not a diagonal matrix

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);
% count = how many zeros not on the diagonal
count=0;
for i=1:1:n
for j=1:1:n
if A(i,j)==0 & i~=j
count=count+1;
end
end
end
disp(‘FIRST WAY’)
if count==n^2-n
disp(‘Matrix is diagonal’)
else
disp(‘Matrix is NOT diagonal’)
end

%% SECOND SOLUTION
% This is based on finding if any of the off-diagonal elements
% are nozero. As soon as this condition is met, the matrix can be
% deemed not diagonal. If the condition is never met, the matrix is
% diagonal

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);
% flag = keeps track if it is diagonal or not
% flag = 1 if matrix is diagonal
% flag = 2 if matrix is not diagonal

% Assuming matrix is diagonal
flag=1;
for i=1:1:n
for j=1:1:n
% flag=2 if off-diagonal element is nonzero.
if A(i,j)~=0 & i~=j
flag=2;
end
end
end
disp(‘ ‘)
disp(‘SECOND WAY’)
if flag==1
disp(‘Matrix is diagonal’)
else
disp(‘Matrix is NOT diagonal’)
end

%% THIRD SOLUTION
% This is based on finding if the sum of the absolute value of
% the off-diagonal elements is nonzero.
% If the sum is nonzero, the matrix is NOT diagonal.
% If the sum is zero, the matrix is diagonal

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);

% sum_off_diagonal= sum of absolute value of off-diagonal elements
sum_off_diagonal=0;
for i=1:1:n
for j=1:1:n
if i~=j
sum_off_diagonal=sum_off_diagonal+abs(A(i,j));
end
end
end

disp(‘ ‘)
disp(‘THIRD WAY’)
if sum_off_diagonal==0
disp(‘Matrix is diagonal’)
else
disp(‘Matrix is NOT diagonal’)
end

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.