## How do I solve a boundary 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 a boundary 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 a boundary value ordinary differential equation?

%% SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://nm.mathforcollege.com/blog/ode_boundaryvalue.m;
% Last Revised : May 26, 2009;
% Abstract: This program shows you how to solve an
%           boundary value ordinary differential equation.
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to solve’)
disp(‘   a boundary value ordinary differential equation’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/ode_boundaryvalue.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   May 26, 2009’)
disp(‘ ‘)

%% INPUTS
% Solve the ordinary differential equation
% y”-0.0000002y=0.000000075*x*(75-x)
% Define x as a symbol
syms x
%The ODE
ode_eqn=’D2y-0.000002*y=0.000000075*x*(75-x)’;
% The initial conditions
iv_1=’y(0)=0′;
iv_2=’y(75)=0′;
% The value at which y is sought at
xval=37.5;
%% DISPLAYING INPUTS

disp(‘INPUTS’)
func=[‘  The ODE to be solved is ‘ ode_eqn];
disp(func)
iv_explain=[‘  The boundary 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://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

## 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://nm.mathforcollege.com/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://nm.mathforcollege.com/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://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

## 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 differentiate 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 differentiate a 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 differentiate a function?

%% SUMMARY
% Language : Matlab 2008a
% Authors : Autar Kaw
% Mfile available at
% http://nm.mathforcollege.com/blog/differentiation.m
% Last Revised : March 21, 2009
% Abstract: This program shows you how to differentiate a given function

%% INTRODUCTION
clc
clear all
disp(‘ABSTRACT’)
disp(‘   This program shows you how to differentiate’)
disp(‘   a given function and then find its value’)
disp(‘   at a given point’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com&#8217;)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/differentiation.m&#8217;)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   March 21, 2009’)
disp(‘ ‘)

%% INPUTS

% Differentiate 7 exp(3*x) once and find the value of the
% first derivative at x=0.5
% Define x as a symbol
syms x
% Defining the function to be differentiated
y=7*exp(3*x);
% Defining the point where you want to find the derivative
xx=0.5;

%% DISPLAYING INPUTS
disp(‘INPUTS’)
func=[‘  The function is to be differentiated is ‘ char(y)];
disp(func)
fprintf(‘  Value of x where you want to find the derivative, x= %g’,xx)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE
% Finding the derivative using the diff command
% Argument 1 is the function to be differentiated
% Argument 2 is the variable with respect to which the
%    function is to be differentiated – the independent variable
% Argument 3 is the order of derivative
dydx=diff(y,x,1);
% subs command substitues the value of x
dydx_val=subs(dydx,x,xx);
%% DISPLAYING OUTPUTS
disp(‘OUTPUTS’)
derivative_func=[‘  The derivative of function ‘ char(y) ‘ is ‘ char(dydx)];
disp(derivative_func)
fprintf(‘  Value of dydx at x=%g is =%g’,xx,dydx_val)
disp(‘  ‘)

______________________________________________________________________

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

## 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.

## A Matlab program for comparing Runge-Kutta methods

In a previous post, we compared the results from various 2nd order Runge-Kutta methods to solve a first order ordinary differential equation. In this post, I am posting the matlab program. It is better to download the program as single quotes in the pasted version do not translate properly when pasted into a mfile editor of MATLAB or see the html version for clarity and sample output .

Do your own testing on a different ODE, a different value of step size, a different initial condition, etc. See the inputs section below that is colored in bold brick.

% Simulation : Comparing the Runge Kutta 2nd order method of
% solving ODEs
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : July 12, 2008
% Abstract: This program compares results from the
% exact solution to 2nd order Runge-Kutta methods
% of Heun’s method, Ralston’s method, Improved Polygon
% method, and directly using the three terms of Taylor series
clc
clear all
clc
clf
disp(‘This program compares results from the’)
disp(‘exact solution to 2nd order Runge-Kutta methods’)
disp(‘of Heuns method, Ralstons method, Improved Polygon’)
disp(‘ method, and directly using the three terms of Taylor series’)

%INPUTS. If you want to experiment these are only things
% you should and can change. Be sure that the ode has an exact
% solution
% Enter the rhs of the ode of form dy/dx=f(x,y)
fcnstr=’sin(5*x)-0.4*y’ ;
% Initial value of x
x0=0;
% Initial value of y
y0=5;
% Final value of y
xf=5.5;
% number of steps to go from x0 to xf.
% This determines step size h=(xf-x0)/n
n=10;

%REST OF PROGRAM
%Converting the input function to that can be used
f=inline(fcnstr) ;

% EXACT SOLUTION
syms x
eqn=[‘Dy=’ fcnstr]
% exact solution of the ode
exact_solution=dsolve(eqn,’y(0)=5′,’x’)
% geting points for plotting the exact solution
xx=x0:(xf-x0)/100:xf;
yy=subs(exact_solution,x,xx);
yexact=subs(exact_solution,x,xf);
plot(xx,yy,’.’)
hold on

% RUNGE-KUTTA METHODS
h=(xf-x0)/n;
% Heun’s method
a1=0.5;
a2=0.5;
p1=1;
q11=1;
xr=zeros(1,n+1);
yr=zeros(1,n+1);
%Initial values of x and y
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=f(xr(i),yr(i));
k2=f(xr(i)+p1*h,yr(i)+q11*k1*h);
yr(i+1)=yr(i)+(a1*k1+a2*k2)*h;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_heun=yr(n+1);
% Absolute relative true error for value using Heun’s Method
et_heun=abs((y_heun-yexact)/yexact)*100;
hold on
xlabel(‘x’)
ylabel(‘y’)
title_name=[‘Comparing exact and Runge-Kutta methods with h=’ num2str(h)] ;
title(title_name)
plot(xr,yr, ‘color’,’magenta’,’LineWidth’,2)
% Midpoint Method (also called Improved Polygon Method)
a1=0;
a2=1;
p1=1/2;
q11=1/2;
%Initial values of x and y
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=f(xr(i),yr(i));
k2=f(xr(i)+p1*h,yr(i)+q11*k1*h);
yr(i+1)=yr(i)+(a1*k1+a2*k2)*h;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_improved=yr(n+1);
% Absolute relative true error for value using Improved Polygon Method
et_improved=abs((y_improved-yexact)/yexact)*100;
hold on
plot(xr,yr,’color’,’red’,’LineWidth’,2)

% Ralston’s method
a1=1/3;
a2=2/3;
p1=3/4;
q11=3/4;
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=f(xr(i),yr(i));
k2=f(xr(i)+p1*h,yr(i)+q11*k1*h);
yr(i+1)=yr(i)+(a1*k1+a2*k2)*h;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_ralston=yr(n+1);
% Absolute relative true error for value using Ralston’s Method
et_ralston=abs((y_ralston-yexact)/yexact)*100;
hold on
plot(xr,yr,’color’,’green’,’LineWidth’,2)

% Using first three terms of the Taylor series
syms x y;
fs=char(fcnstr);
% fsp=calculating f'(x,y) using chain rule
fsp=diff(fs,x)+diff(fs,y)*fs;
%Initial values of x and y
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=subs(fs,{x,y},{xr(i),yr(i)});
kk1=subs(fsp,{x,y},{xr(i),yr(i)});
yr(i+1)=yr(i)+k1*h+1/2*kk1*h^2;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_taylor=yr(n+1);
% Absolute relative true error for value using Taylor series
et_taylor=abs((y_taylor-yexact)/yexact)*100;
hold on
plot(xr,yr,’color’,’black’,’LineWidth’,2)
hold off
legend(‘exact’,’heun’,’midpoint’,’ralston’,’taylor’,1)

% THE OUTPUT
fprintf(‘\nAt x = %g ‘,xf)
disp(‘ ‘)
disp(‘_________________________________________________________________’)
disp(‘Method Value Absolute Relative True Error’)
disp(‘_________________________________________________________________’)
fprintf(‘\nExact Solution %g’,yexact)
fprintf(‘\nHeuns Method %g %g ‘,y_heun,et_heun)
fprintf(‘\nImproved method %g %g ‘,y_improved,et_improved)
fprintf(‘\nRalston method %g %g ‘,y_ralston,et_ralston)
fprintf(‘\nTaylor method %g %g ‘,y_taylor,et_taylor)
disp( ‘ ‘)

________________________________________________________________________________________________

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