Using MATLAB To Solve Differential Equations

Transcription

ECE 350 – Linear Systems IMATLAB Tutorial #3Using MATLAB to Solve Differential EquationsThis tutorial describes the use of MATLAB to solve differential equations. Two methodsare described. The first uses one of the differential equation solvers that can be calledfrom the command line. The second uses Simulink to model and solve a differentialequation.Solving First Order Differential Equations with ode45The MATLAB commands ode 23 and ode 45 are functions for the numerical solution ofordinary differential equations. They use the Runge-Kutta method for the solution ofdifferential equations. This example uses ode45. The function ode45 uses higher orderformulas and provides a more accurate solution than ode 23.In this example we will solve the first order differential equation:dy 2y u(t) u(t 1)dtover the range 0 t 5. We will assume that the initial value y(0) 0. The first step isto write the equation in the form:dy 2y u(t) u(t 1)dtNext we need to create a MATLAB m-file to calculate dy/dt. From the MATLAB Homemenu select: New ! Function. The basic structure of a function file will appear in theeditor window. Modify the file to read:function dy f(t,y)dy -2*y (t 0)-(t 1);This creates a function that calculates the value of dy (actually dy/dt). The function iscalled “f” and the file should be named: f.m when it is saved since all MATLABfunctions should be saved in a file with the same name as the function. Convince

yourself that the right hand side of the equation for dy corresponds to the expressionabove. Make sure that you save this function in your working directory or one that is inyour path.Once the function for dy/dt has been created we need to specify the end points of the timerange (0 t 5) and the initial value of y (zero). This is done in the MATLAB commandwindow with the following commands: t [0 5]; inity 0;The ode45 command can now be called to compute and return the solution for y alongwith the corresponding values of t using: [t,y] ode45(@f, t, inity);Note that the @f refers to the name of the function containing the function. The solutioncan be viewing with: plot(t,y)as shown in figure 1.Figure 1. Solution of First Order Differential EquationSolving First Order Differential Equations with ode45The ode45 command solves first order differential equations. In order to use thiscommand to solve a higher order differential equation we must convert the higher orderequation to a system of first order differential equations. For example, in order to solvethe second order equation:

d2 ydy 3 2y cos 2tdt 2dtwith y(0) 0 anddy 1dt t 0for 0 t 10 we must define a variable x such that:x dydtThus we can rewrite the differential equation as:d2 ydy 3 2y cos 2t2dtdtor in terms of x:dx 3x 2y cos 2tdtWe will use the ode45 command to solve the system of first order differential equations:dx 3x 2y cos 2tdtdy xdtWe place these two equations into a column vector, z, where:! x z #&#" y &%and!#dz # dt ##"dxdtdydt &&&&%A function is created in a MATLAB m-file to define dz/dt named f2 which contains:function dz f2(t,z)dz [-3*z(1)-2*z(2) cos(2*t); z(1)];Note that the first row of dz contains the equation for dx/dt and the second an equationfor dy/dt. The values of x and y are referred to as z(1) and z(2) since they are the firstand second elements of the vector z. The range of t values and the initial conditions areinput at the command line with: t [0 10]; initz [1;0];

and the equations are solved and plotted using: [t,z] ode45(@f2, t, initz); plot(t, z(:,2))Figure 2. Solution of Second Order Differential EquationModeling Linear Systems Using SimulinkSimulink is a companion program to MATLAB and is included with the student version.It is an interactive system for simulating linear and nonlinear dynamic systems. It is agraphical mouse-driven program that allows you to model a system by drawing a blockdiagram on the screen and manipulating it dynamically. It can work with linear,nonlinear, continuous time, discrete time, multivariable, and multirate systems. In thissection the basic operation of Simulink will be described. The next section willdemonstrate the use of Simulink to solve differential equations.1. Open MATLAB and in the command window, type: simulink at the prompt.2. After a few seconds Simulink will open and the Simulink Library Browser will openas shown in figure 3. It is important to note that the list of libraries may be differenton your computer. The libraries are a function of the toolboxes that you haveinstalled.

Figure 3. Simulink Library Browser3. Click on the Create a New Model icon in the Library Browser window. An additionalwindow will open. This is where you will build your Simulink models.4. Click on the arrow next to “Simulink” in the Library Browser. A list of sub-librarieswill appear including Continuous, Discrete, etc. These sub-libraries contain the mostcommon Simulink blocks.5. Click once on the “Sources” sub-library. You should see a listing of blocks as shownin the right column of figure 4.

Figure 4. Source Blocks in the Simulink Library6. Scroll down this list until you see the Sine Wave icon. Click once on the icon to selectit and drag this icon into the blank model window.7. Click once on the “Sinks” sub-library in the left part of the Library Browser. Clickand drag the “Scope” icon to the model window to the right of the Sine Wave block.The model window should now appear as shown in figure 5. Make sure you haveused “Scope”, not “Floating Scope”.

Figure 5. Model Window with Sine Wave and Scope Blocks8. Next we want to connect the Sine Wave to the Scope block. Move the mouse overthe output terminal of the Sine Wave block until it becomes a crosshair. Click anddrag the wire to the input terminal of the Scope block. The cursor will become adouble cursor when it is in the correct position. Release the mouse button and thewire will snap into place. Your completed system should now appear as shown infigure 6.Figure 6. Completed System for Viewing Sine Wave

9. In the model window click on the green arrow “Run” icon. You will hear a beepwhen the simulation is complete. Double click on the Scope icon and it will open anddisplay the output of the sine wave block. Try clicking on the various icons in theScope window. You should be able to cause the axes to readjust to display thewaveform as shown in figure 7.Figure 7. Scope Display10. Note that the period of the sine wave is just over six. What frequency does thiscorrespond to? Return to the model window and double click on the Sine Waveblock. The Sine Wave parameters window will open as shown in figure 8.

Figure 8. Sine Wave Block Parameters11. The frequency was set to 1 rad/sec. Any of the parameters shown can be altered tochange the sine wave. Change the entry in the frequency parameter to: 2*pi*10.This is a 10 Hz sine wave. Click OK and re-simulate. Use the icons to expand thescope display. What do you observe? The sine wave should NOT be displayedproperly. This is because at this higher frequency the sampling rate has not been sethigh enough.12. Return to the model window and select the Model Configuration Parameters icon.The Stop Time for the simulation defaults to 10. Change it to 1, click OK, and re-runthe simulation. You should now be able to properly view the sine wave on the scope.Verify that the period of each cycle is one 0.1 second as you expect for a 10 Hz wave.13. Another method for modifying the display is to specify the step size. Return to themodel window and select the Configuration Parameters window again. Note that theSolver Type defaults to “variable-step”. Change the Solver Type to “fixed step” andthen enter the value 0.01 into the Fixed Step Size entry. Click OK and run thesimulation once again. Note that the display is now smooth.

14. Drag the Pulse Generator from the Source sub-library into the model window.Double click the Pulse Generator block and modify the parameters as shown in figure9 below. Click OK.Figure 9. Pulse Generator Parameters15. In the scope window, click on the Parameters icon (second from the left). Change thenumber of axes to 2 and click OK. Return to the model window and note that thescope now has two inputs. Connect a wire from the Pulse Generator to the secondinput. Your model should now appear as shown in figure 10.

Figure 10. Model Window with Two Sources16. Re-run the simulation and observe the two waveforms as shown in figure 11 below.Figure 11. Two Input Scope Window17. Open the Continuous sub-library. Drag the Integrator block into the model windowand reconfigure the diagram as shown in figure 12.

Figure 12. Simulink Diagram18. Re-simulate the system and observe the output on the scope. Note that the secondtrace is now the integral of the pulse train input.Modeling a Differential Equation in SimulinkNow we will use Simulink to model the differential equations solved earlier. Beginningwith the first order differential equation:dy 2y u(t) u(t 1)dtFirst we need to simulate the single one-second pulse. Since the simulation will run for0 t 5, we can use the pulse generator to generate a five second pulse with a 1 secondpulse width. Since it will not repeat within the five seconds, this simulates a single pulse.Remove the integrator from the model shown in figure 12 and connect the pulsegenerator directly to the scope input. Change the pulse generator to have the settings asshown in figure 13. The 20% pulse width corresponds to one second as needed.Change the simulation time in the configuration parameters to five seconds and simulatethe system. Specify fixed-step samples of 0.01 seconds. Verify that the pulse generatoris outputting the desired pulse for the differential equation above.

Figure 13. Pulse Generator Settings for 1 Second PulseThe differential equation above can be written as:dy 2y u(t) u(t 1) 2y p(t)dtwhere p(t) is the one second pulse. The right hand side of this equation can be modeledin Simulink using the model shown in figure 14. The subtraction block and the gainblock are found in the Math Operations sub-library.Figure 14. Simulink Model

The labels on the wires are inserted by double clicking on the wires and typing in the text.By double clicking on the gain block you can set its constant to two. If the input to thegain block is y, then the output of the subtractor is dy/dt. By passing this output throughan integrator, the input y is found. The model in figure 15 should now simulate theoriginal differential equation.Figure 15. Simulink Model of First Order Differential EquationRun the simulation for 5 seconds and confirm that the scope displays the same output aswe obtained in the first section of this tutorial.We can also simulate the second order differential equation:d2 ydy 3 2y cos 2t2dtdtin Simulink. This equation can be re-written as:d2 ydy 3 2y cos 2t2dtdtand modeled in Simulink as shown in figure 16. In order to get the three input subtractor,use the two input subtractor selected above. Double click on the block and change the

“List of Signs” to: - - . This will result in a block that adds the top input and subtractsthe second and third inputs.Figure 16. Simulink ModelIn order to get y and dy/dt we need to integrate the output of this summer twice. If weadd two integrators to the output of figure 16 we can feed the values of y and dy/dt backto the inputs as shown in figure 17.Figure 17. Simulink Model of 2nd Order Differential EquationSimulate the circuit for 10 seconds. The output shown in figure 18 is obtained on thescope.

Figure 18. MATLAB Scope OutputIf we want to compare this output to that obtained in the earlier solution we can save theoutput to the workspace. The simout block found in the Sinks sub-library is added asshown in figure 19. Double click the simout block and change the format of the output toan Array.Figure 19. Simulink Model of 2nd Order CircuitSimulate the model again. Note that an array named simout appears in the workspace.Plot this output along with the output obtained for the command line solution of thesecond order differential equation. The following commands are used: z [1; 0]; t [0 10]; [t, z] ode45(@f2, t, z);

plot(t,z(:,2)) hold on plot(0:.01:10, simout, 'r')The two solutions appear in figure 20. Note that the time values for the simout plot aredetermined by the fixed spacing of 0.01 specified in the configuration parameters inSimulink.Figure 20. Two Solutions of the 2nd Order Differential EquationThe solutions are not the same. There are two reasons. First of all, the sinusoidal input isa “sine” and not a “cos”. In order to convert a sine to a cos we need to apply a phase shiftof π/2. Double click on the Sine Wave block and enter: pi/2 for the phase. Re-simulateand plot along with the command line solution as shown in figure 21.Figure 21. Two Solutions of the 2nd Order Differential Equation

Now the steady state portions of the solutions are the same. However, the initial transientresponse is not. This is because the two solutions were obtained with different initialconditions. In the earlier command line solution the initial conditions: y(0) 0 anddy 1 were used. The integrator blocks in Simulink allow for the inclusion ofdt t 0initial conditions. The second integrator outputs the value of y. Thus, the default initialcondition of zero is correct. However the first integrator outputs dy/dt. Double click onthe first integrator and change the initial condition to one. After simulating and replotting both solutions the curves in figure 22 are obtained. Note that the solutions nowmatch. This confirms that either method will provide a correct solution.Figure 22. Two Solutions of the 2nd Order Differential Equation

Next we need to create a MATLAB m-file to calculate dy/dt. From the MATLAB Home menu select: New ! Function. The basic structure of a function file will appear in the editor window. Modify the file to read: function dy f(t,y) dy -2*y (t 0)-(t 1); This creates a function that calculates the value of dy (actually dy/dt). The function is