R TUTORIAL for the First Course. Part 3: Numerical Methods

Prof. Vladimir A. Dobrushkin

This tutorial contains many R scripts. You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to

Email Vladimir Dobrushkin

Introduction
Solving ODEs Numerically
We present a couple of examples to solve an initial value problem for the first order differential equation in normal form using standard Matlab subroutine. 

function solveDifferentialEquation

              %This function solves y'+y = x from starttime to endtime with y(0) = intialCondition


%Parameter initializations
              starttime = 0;
              endtime = 20;
              initialCondition = -2;
 function yprime = DifferentialEquationFunction(x,y)
              yprime = y - x;
              end
[xVals,yVals] = ode45(@DifferentialEquationFunction,[starttime,endtime],initialCondition);
              plot(xVals,yVals);
              end %Function Complete


			function solveDifferentialEquation2
              % Function to solve the differential equation
              % of a particle in motion due to impact of an external sinusoidal
              % force ( yprime = F*sin(w*t) - C*y )
              % from a starttime to stoptime, where the intial y value is assigned as
              % initialCondition

              %Parameter initializations
              n = 5;
              onesv = [1, 1.45, 1.34, 1.81, 3.5]; %acts as a scalar vector that scales the parameters by its %entries.
              C = 10.*onesv;
              F = 1.*onesv;
              w = 1.*onesv; %With scaled frequencies,w, the plots should show a trend of increasing spread with t5he 3rd plot as an exception since 1.34 defies the pattern.



starttime = 0;
              endtime = 20;
              initialCondition = 0;
              counter = 1;

          %Loop to solve the ODE and plot five different graphs.

              for indeX = 1:n
              [t,y] = ode45(@DifferentialEquationFunction2,[starttime,endtime],initialCondition);
              counter = counter + 1;
              figure(indeX)
              plot(t,y);

              title('Solution to differential equation')
              xlabel('time')
              ylabel('y-value')
              end
            
function yprime = DifferentialEquationFunction2(t,y)      % Function defining the ODE
              yprime = F(counter)*sin(w(counter)*t) - C(counter)*y;
              end

              end          %Function Complete

There is a Matlab subroutine, credited to John Polking, called odesolve.m

Example 1.1.1:

Example 1.1.2:

II. Plotting