{"id":224,"date":"2010-03-25T18:21:31","date_gmt":"2010-03-25T18:21:31","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2010\/03\/25\/solving-ordinary-differential-equations\/"},"modified":"2020-03-26T11:56:38","modified_gmt":"2020-03-26T16:56:38","slug":"solving-ordinary-differential-equations","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2010\/03\/25\/solving-ordinary-differential-equations\/","title":{"rendered":"Solving Ordinary Differential Equations"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>I have recently handled several help requests for solving differential equations in MATLAB.  All of the cases I worked on\r\n         boil down to how to transform the higher-order equation(s) given to a system of first order equations.  In this post I will\r\n         outline how to accomplish this task and solve the equations in question.  What I am <i>not<\/i> going to talk about is details of ODE solving algorithms (maybe another time).\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">Simple System<\/a><\/li>\r\n         <li><a href=\"#2\">Conditions, Initial and Otherwise<\/a><\/li>\r\n         <li><a href=\"#3\">Transform Equation<\/a><\/li>\r\n         <li><a href=\"#4\">Try It!<\/a><\/li>\r\n         <li><a href=\"#11\">Transforming Equations<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Simple System<a name=\"1\"><\/a><\/h3>\r\n   <p>Let's start with a simple, linear, second order system, a mass-spring system with motion constrained to one-dimension, a horizontal\r\n      one (and therefore no gravity).  If <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq81831.png\">  is the mass and <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq86607.png\">  is the spring constant, the equations of motion for the system are:\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq14881.png\"> <\/p>\r\n   <h3>Conditions, Initial and Otherwise<a name=\"2\"><\/a><\/h3>\r\n   <p>To solve this system, we need to know <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq86924.png\"> , <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq81831.png\"> , <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq86607.png\"> , and initial conditions, e.g., <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq89740.png\">  (also known as position and velocity). Let's simplify things and set <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq78918.png\"> , i.e., no external forces.  Let's also set some initial conditions, <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq66156.png\"> , in other words, start with the spring unstretched and the mass moving. I end up with this system:\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq07197.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq66156.png\"> <\/p>\r\n   <h3>Transform Equation<a name=\"3\"><\/a><\/h3>\r\n   <p>Looking in the help, I need to set up an system of equations to enable me to use one of the <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2010a\/techdoc\/ref\/ode23tb.html\">numerical ODE solvers<\/a> in MATLAB.\r\n   <\/p>\r\n   <p>To start the transformation, let me define a new variable that I will substitute in the system.<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq09728.png\"> <\/p>\r\n   <p>I can derive <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq08970.png\"> <\/p>\r\n   <p>and now rewrite my ODE system in terms of <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq66157.png\"> .\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq09728.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq89615.png\"> <\/p>\r\n   <p>with initial conditions as above.<\/p>\r\n   <p>Now let me reorganize these 2 equations in a vector\/matrix equation where <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq67049.png\"> <\/p>\r\n   <p>I now write me equation solely in terms of <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq88768.png\"> , the new vector (consisting of position and velocity).  With that in mind, I will reorganize the existing equations first\r\n      so I have <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq02286.png\">  on the left-hand sides.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq73917.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq75380.png\"> <\/p>\r\n   <p>or, in terms of <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq88768.png\"> ,\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq70476.png\"> <\/p>\r\n   <h3>Try It!<a name=\"4\"><\/a><\/h3>\r\n   <p>Let's try it. Set values for <tt>m<\/tt> and <tt>k<\/tt><\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">m = 1;\r\nk = 10;<\/pre><p>I can create the ODE code in a file, or I can set up the equations as an anonymous function which is what I'll do here.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">springmass = @(t,z)[z(2); -k * z(1)\/ m];<\/pre><p>Set up the initial conditions.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">ic = [0; 1];<\/pre><p>Solve between <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_eq69326.png\"> .\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">tspan = [0 10];<\/pre><p>Call an ODE solver and plot the results.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">[t,y] = ode23(springmass, tspan, ic);\r\nplot(t,y(:,1)), title(<span style=\"color: #A020F0\">'Position vs. Time'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_01.png\"> <p>And the velocity.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">plot(t,y(:,2)), title(<span style=\"color: #A020F0\">'Velocity vs. Time'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/224\/odeSolve_02.png\"> <p>You can see the left-most points on the plots match the initial conditions specified.<\/p>\r\n   <h3>Transforming Equations<a name=\"11\"><\/a><\/h3>\r\n   <p>Do you need to recast problems to fit into a particular formulation such as for solving ODEs?  Let's hear about them <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=224#respond\">here<\/a>.\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_ec707d92295e44a7add99c2e78318fbc() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='ec707d92295e44a7add99c2e78318fbc ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' ec707d92295e44a7add99c2e78318fbc';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Loren Shure';\r\n        copyright = 'Copyright 2010 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_ec707d92295e44a7add99c2e78318fbc()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.10<br><\/p>\r\n<\/div>\r\n<!--\r\nec707d92295e44a7add99c2e78318fbc ##### SOURCE BEGIN #####\r\n%% Solving Ordinary Differential Equations\r\n% I have recently handled several help requests for solving differential\r\n% equations in MATLAB.  All of the cases I worked on boil down to how to\r\n% transform the higher-order equation(s) given to a system of first order\r\n% equations.  In this post I will outline how to accomplish this task and\r\n% solve the equations in question.  What I am _not_ going to talk about is\r\n% details of ODE solving algorithms (maybe another time).\r\n%% Simple System\r\n% Let's start with a simple, linear, second order system, a mass-spring\r\n% system with motion constrained to one-dimension, a horizontal one (and\r\n% therefore no gravity).  If $m$ is the mass and $k$ is the spring\r\n% constant, the equations of motion for the system are:\r\n%\r\n% $F = m * x(t)'' + k * x(t)$\r\n%\r\n%% Conditions, Initial and Otherwise\r\n% To solve this system, we need to know $F$, $m$, $k$, and initial\r\n% conditions, e.g., $x(0), x'(0)$ (also known as position and velocity).\r\n% Let's simplify things and set $F=0$, i.e., no external forces.  Let's\r\n% also set some initial conditions, $x(0) = 0, x'(0) = 1$, in other words,\r\n% start with the spring unstretched and the mass moving. I end up with this\r\n% system:\r\n%\r\n% $m * x(t)'' + k * x(t) = 0$\r\n%\r\n% $x(0) = 0, x'(0) = 1$\r\n%\r\n%% Transform Equation\r\n% Looking in the help, I need to set up an system of equations to enable me\r\n% to use one of the\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2010a\/techdoc\/ref\/ode23tb.html\r\n% numerical ODE solvers> in MATLAB.\r\n%\r\n% To start the transformation, let me define a new variable that I will\r\n% substitute in the system.\r\n%\r\n% $y = x'$\r\n%\r\n% I can derive\r\n% $y' = x''$\r\n%\r\n% and now rewrite my ODE system in terms of $x, y, y'$.\r\n%\r\n% $y = x'$\r\n%\r\n% $m * y(t)' + k * x(t) = 0$\r\n%\r\n% with initial conditions as above.\r\n%\r\n% Now let me reorganize these 2 equations in a vector\/matrix equation where\r\n% $z = [ x ; y ]$\r\n%\r\n% I now write me equation solely in terms of $z$, the new vector\r\n% (consisting of position and velocity).  With that in mind, I will\r\n% reorganize the existing equations first so I have $x', y'$ on the\r\n% left-hand sides.\r\n%\r\n% $x' = y$\r\n%\r\n% $y' =  -k * x(t) \/ m$\r\n%\r\n% or, in terms of $z$,\r\n%\r\n% $z' = [z(2); -k * z(1)\/ m]$\r\n%% Try It!\r\n% Let's try it. Set values for |m| and |k|\r\nm = 1;\r\nk = 10;\r\n%%\r\n% I can create the ODE code in a file, or I can set up the equations\r\n% as an anonymous function which is what I'll do here.\r\nspringmass = @(t,z)[z(2); -k * z(1)\/ m];\r\n%%\r\n% Set up the initial conditions.\r\nic = [0; 1];\r\n%%\r\n% Solve between $t = 0, t = 10$.\r\ntspan = [0 10];\r\n%%\r\n% Call an ODE solver and plot the results.\r\n[t,y] = ode23(springmass, tspan, ic);\r\nplot(t,y(:,1)), title('Position vs. Time')\r\n%%\r\n% And the velocity.\r\nplot(t,y(:,2)), title('Velocity vs. Time')\r\n%%\r\n% You can see the left-most points on the plots match the initial\r\n% conditions specified.\r\n%% Transforming Equations\r\n% Do you need to recast problems to fit into a particular formulation such\r\n% as for solving ODEs?  Let's hear about them\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=224#respond here>.\r\n\r\n\r\n\r\n##### SOURCE END ##### ec707d92295e44a7add99c2e78318fbc\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      I have recently handled several help requests for solving differential equations in MATLAB.  All of the cases I worked on\r\n         boil down to how to transform the higher-order... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2010\/03\/25\/solving-ordinary-differential-equations\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[87,39],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/224"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=224"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/224\/revisions"}],"predecessor-version":[{"id":3619,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/224\/revisions\/3619"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=224"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=224"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=224"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}