Solving Set of Second Order ODEs with Matlab ODE45 function

Question

Introduction

NOTE IN CODE AND DISUSSION:

A single d is first derivative A double d is second derivative

I am using Matlab to simulate some dynamic systems through numerically solving the governing LaGrange Equations. Basically a set of Second Order Ordinary Differential Equations. I am using ODE45. I found a great tutorial from Mathworks (link for tutorial below) on how to solve a basic set of second order ordinary differential equations.

https://www.mathworks.com/academia/student_center/tutorials/source/computational-math/solving-ordinary-diff-equations/player.html

Based on the tutorial I simulated the motion for an elastic spring pendulum by obtaining two second order ordinary differential equations (one for angle theta and the other for spring elongation) shown below:

theta double prime equation:

M*thetadd*(L + del)^2 + M*g*sin(theta)*(L + del) + M*deld*thetad*(2*L + 2*del) = 0

del (spring elongation) double prime equation:

K*del + M*deldd - (M*thetad^2*(2*L + 2*del))/2 - M*g*cos(theta) = 0

Both equations above have form ydd = f(x, xd, y, yd)

I solved the set of equations by a common reduction of order method; setting column vector z to [theta, thetad, del, deld] and therefore zd = [thetad, thetadd, deld, deldd]. Next I used two matlab files; a simulation file and a function handle file for ode45. See code below of simulation file and function handle file:

Simulation File

%ElasticPdlmSymMainSim
clc
clear all;
%Define parameters
global M K L g;
M = 1;
K = 25.6;
L = 1;
g = 9.8;
% define initial values for theta, thetad, del, deld
theta_0 = 0;
thetad_0 = .5;
del_0 = 1;
deld_0 = 0;
initialValues = [theta_0, thetad_0, del_0, deld_0];
% Set a timespan
t_initial = 0;
t_final = 36;
dt = .01;
N = (t_final - t_initial)/dt;
timeSpan = linspace(t_final, t_initial, N);
% Run ode45 to get z (theta, thetad, del, deld)
[t, z] = ode45(@OdeFunHndlSpngPdlmSym, timeSpan, initialValues);

Here is the function handle file:

function dz = OdeFunHndlSpngPdlmSym(~, z)
% Define Global Parameters
global M K L g
% Take output from SymDevFElSpringPdlm.m file for fy1 and fy2 and
% substitute into z2 and z4 respectively
% z1 and z3 are simply z2 and z4 
% fy1=thetadd=z(2)= -(M*g*sin(z1)*(L + z3) + M*z2*z4*(2*L + 2*z3))/(M*(L + z3)^2)
% fy2=deldd=z(4)=((M*(2*L + 2*z3)*z2^2)/2 - K*z3 + M*g*cos(z1))/M
% return column vector [thetad; thetadd; deld; deldd]
dz = [z(2);
    -(M*g*sin(z(1))*(L + z(3)) + M*z(2)*z(4)*(2*L + 2*z(3)))/(M*(L + z(3))^2);
    z(4);
    ((M*(2*L + 2*z(3))*z(2)^2)/2 - K*z(3) + M*g*cos(z(1)))/M];

Question

However, I am coming across systems of equations where the variables can not be solved for explicitly as is the case with spring pendulum example. For one case I have the following set of ordinary differential equations:

y double prime equation

ydd - .5*L*(xdd*sin(x) + xd^2*cos(x) + (k/m)*y - g = 0

x double prime equation

 .33*L^2*xdd - .5*L*ydd*sin(x) - .33*L^2*C*cos(x) + .5*g*L*sin(x) = 0

L, g, m, k, and C are given parameters.

Note that x'' term appears in y'' equation and y'' term appears in x'' equation so I am not able to use reduction of order method. Can I use Matlab ODE45 to solve the set of ordinary differential equations in the second example in a manner similar to first example?

Thanks!

You're missing a parenthesis somewhere in your "y double prime equation". — horchler
That equation is for question purposes and not part of any code. — PatStarks
Duplicate or follow-up to stackoverflow.com/questions/41825370/… — Lutz Lehmann
Well this is a clarification because in that code I did not mention that x and y variables are functions of time and that the goal is to obtain their values as time marches forward based on initial conditions. — PatStarks

Jonathan Chiang Jonathan Chiang · Accepted Answer · 2017-01-25T19:18:42

This problem can be solved by working out some of the math by hand. The equations are linear in xdd and ydd so it should be straightforward to solve.

ydd - .5*L*(xdd*sin(x) + xd^2*cos(x)) + (k/m)*y - g = 0

.33*L^2*xdd - .5*L*ydd*sin(x) - .33*L^2*C*cos(x) + .5*g*L*sin(x) = 0

can be rewritten as

-.5*L*sin(x)*xdd +             ydd = -.5*L*xd^2*cos(x) - (k/m)*y + g

     .33*L^2*xdd - .5*L*sin(x)*ydd = .33*L^2*C*cos(x) - .5*g*L*sin(x)

which is the form A*x=b.

For more complex systems, you can look into the fsolve function.

Solving Set of Second Order ODEs with Matlab ODE45 function

1 Answers