lab 12

pdf

School

San Jose State University *

*We aren’t endorsed by this school

Course

130

Subject

Mechanical Engineering

Date

Jan 9, 2024

Type

pdf

Pages

13

Uploaded by PrivateRook16679

Report
Lab 12: Using MATLAB's Solvers for ODEs © Matthew Leineweber, PhD BME 130 - Numerical Methods in Biomedical Engineering San Jose State University Instructor:Prof. Abdulmelik Mohammed Lab TA: Ammar Babiker & Nihar Prakash By: Philippe Moffet, Lab Group #: 03, Year: Fall 2023 Table of Contents Introduction Learning Objectives Ordinary Differential Equations Solving ODEs with Numerical Methods Activity 1 - Visualizing_an ODE Visualizing_the Differential Equation Activity 2 - Plotting_Slope Fields in MATLAB Question Set 1: Activity 3: Plotting_a Slope Field MATLAB's Non-Stiff ODE Solvers: ode23 and ode45 Activity 4: Solving with ode23 and ode45 Activity 5: User-Define Step Size. Solving_Systems of ODEs Activity 6: Solving a System of ODEs Solving Second-Order ODEs Activity 7 Post Lab References and Documentation Introduction Differential equations are essential for solving science and engineering problems, since they are used to model nearly every existing system. While some differential equations can be solved analytically, there are many that require the use of numerical methods. However, there are different numerical methods that are useful for different types of equations, since there is not a “one size fits all” method. Today we will go over solving first-order and second-order ordinary differential equations. Learning Objectives » Understand the steps for solving a first-order ODE » Use canonical form to represent systems of first-order ODEs » Understand how canonical form can be used to solver higher-order ODEs » Learn the difference between MATLAB ODE solvers » Solve first-order ODE’s using MATLAB when given a solution interval and initial conditions Ordinary Differential Equations Recall that ordinary differential equations (ODEs) are defined as mathematical expressions that contains one or more functions of an independent variable (e.g. independent variable x and function y(x)) and its derivatives (y'(x), y"(x), etc.). The most general form for an ODE can be written as: F(x,y,y.y" ..... y") = R(x) However, in most engineering applications we will be most concerned with first and second-order ODEs, where the highest derivatives in the equation are y' and y", respectively. For first order ODEs, the highest derivative in the equation is the first derivative: ir = 0 fx)= E)’(X)
and the general form is %=F@» Similarly, for second-order ODEs, the general form is: d*y , —=F(x,y,y) 12 Y.y So what does it mean to solve an ODE? Essentially, we want to find the function y(x) that satisfies our given differential equation. In many cases, there may be more than one function that meets the requirements of our ODE. In these cases, we need to use initial conditions to narrow our solution down to a single y(x) expression. Solving ODEs with Numerical Methods In Module 4, we explored how we can use the Taylor Series and other approximations to develop “marching” procedures for solving first-order ODEs. Chief among these techniques are the Runge-Kutta methods, which use “increment functions” to represent the local slope of a function over a desired interval. The purpose of this lab is to show you how we can set up and solve single and systems of ODEs using MATLAB’s built-in solvers, many of which rely on Runge-Kutta methods. Setting up ODEs in MATLAB Recall from lecture that numerically solving ODEs requires three key information components: » The expression: y = F(x, y) An interval over which to solve for y: [xq, x| o [nitial conditions: (xq. y(xo)] Therefore, when setting up our ODE in MATLAB, we need to be sure to define each of these components. To illustrate this process, let’s consider an example. Activity 1 - Visualizing an ODE Let’s set up the MATLAB components required to solve the differential equation: y=x+y;. y0)=3 0<x<4 To define the initial condition, we just need to define a variable to represent y(0) = yo =5 Defining our range is slightly more complicated, only in that we cannot represent the continuous range 0 < x < 4 exactly in MATLAB. We need to discretize it, meaning we have to pick a finite number of points over that range for which to solve our equation. We do this with ordered arrays created by either (1) specifying the step size between successive points (i.e. use the colon operator), or (2) defining the number of points desired over the interval (i.e. using linspace). Using the latter approach, we could define: X = linspace(0,4,20); Now that our initial conditions and range are set, we need to define an expression to represent our ODE. We can do this by defining an anonymous function that receives values for x and y as an input, and whose output is the resulting value for y'(x): dydx = @(x,y) x +y; Here we see that the variable dydx contains the function handle that operates on input arguments x and y. That is, dydx will compute the derivative % = y'(x) = x+ y for any (x, y) combination. That’s it! We have now set up all the information we need to begin solving our ODE. Now quickly use our dydx function to evaluate the derivative y'(x) for points (0,0.5) and (1, 2). Report your values to the command window. % <Enter your solution in the space below> yprimel=dydx(0,0.5) yprimel = 0.5000 yprime2=dydx(1,2) yprime2 = 3
Visualizing the Differential Equation Our first order ODE in Activity 1 essentially tells us the slope of the solution function y(x) at any (x, y) point. As such, first-order ODEs are typically plotted as “slope fields”. Using our function for % in MATLAB, we can easily plot the “slope field” (or “direction field”) to get a visual representation of what or solution may look like. For example, the slope field for our ODE in Activity 1 is shown below: Slope Field for = x +y i [ | B | 1 \ \ \ i/ I I I i 1 4 \ \ LRARTE RS { i ! I 1 I 1 | \ \ i/ i ! I I L | i \ \ VAV NASNSSNMNSS YNNI S ¥ S I NNYRNSN VOV VNSNS NSNS R R EREERERERTI. 2N OV OVONN NSNS NSNS \ R T R o rForresrseds s revrvvevevsvhR (1 N N e e e = - & &yeyeyeryreyvevvey o o oey s 7»oi79noua rervrevevvcevveeW o rrrevrcevsve W rrorves s e P NIV Froreere i Fror e el voese syl PR R EEEET PERCEETFF T Figure 1: Slope field representing the solutions to the ordinary differential equation % = x + y. Each (x, y) coordinate contains the derivative of y(x) at that point. Now using these two vectors, we can use MATLAB’s meshgrid function to turn these vectors into 2D arrays, called X and v, that contain all the coordinates of our 3 x 4 grid. Each arrow in the figure denotes the slope of the solution y(x) at that (x, y) point. Visualizing slope fields can be very useful for predicting solutions and confirming whether any solution you calculate fits the original ODE. Let’s explore how we can use MATLAB to plot the slope field. Activity 2 - Plotting Slope Fields in MATLAB Recall that a function y(x) can be visualized as a one-dimensional plot. In other words, there is independent variable (y), which we plotted against one independent variable (x). However, our differential equation depends on two variables: x and y(x). Since y'(x) depends on two variables, and not just one, we will need to use a different plotting command. Refer to Module 2 and Module 3 for examples of plotting 2D and 3D plots in MATLAB. Here, rather than creating surface plots or contour maps, we will use the command quiver to illustrate how the value of y'(x) changes for each (x, y) combination. quiver(x,y,u,v); quiver, as the name suggests, plots arrows (i.e. “vectors”) with horizontal width » and vertical height v at every (x, y) location specified by the input arrays x and y. To use quiver, we need to first define a grid of coordinates. For example, suppose we wanted to create a set of paired (x, y) coordinates corresponding to x = [0,1,2] and y = [0, 1,2,3]. We could represent these points as a 3 x 4 grid. To do so, we would first independently define our and coordinates: 2; 3. J X = 0: y = 0: Now using these two vectors, we can use MATLAB’s meshgrid function to turn these vectors into 2D arrays, called X and Y, that contain all the coordinates of our 3 x 4 grid. [X,Y] = meshgrid(x,y) X = 4x3 (%] 1 2 (%) 1 2 (%] 1 2 (%] 1 2 Y = 4x3
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
w N =R OO w N RO w N P You should see this command produces the following values for X and Y 01 2 00 0 01 2 111 X=lo 12| Y712 2 2 0 1 2 3 3 3 Note that each column in X denotes its horizontal location x, while while each row in Y denotes its vertical location y, so that [X(m,1i),Y(j,n)] represents the point [x;, y;]. Question Set 1: 1. How might we use these new arrays for plotting or defining new functions f(x, y)? We can use meshgrid and dydx = at(x,y)x+y to find v and then use quiver(x,y,u,v) to plot new functions f(x,y) Before creating plot slope field , we need to define the grids of coordinates . Then the two vectors turns to the 2D arrays , they contain all the coordinates of our grid . Those new arrays acts as a specified ( x , y ) location for each plot arrow. 1. Based on the definitions of X and Y, what are the indices corresponding to the point (1,2)? The indices corresponding to the point (1,2)are X=2,Y =3 You can use the find function to figure out where x=1 and y=2 Let’s return now to creating our slope field. We can use meshgrid to create a 2D grid of points to represent our (x, y) points over the range of interest. For visualization, let's recreate the range of shown in Figure 1. Activity 3: Plotting a Slope Field Looking at Figure 1, we can see that it represents the ODE for the range —2 < x < 1 and —-2.5 < y < 2.5. Therefore, we can begin by using meshgrid to define the paired (x, y) coordinates over this range. x = linspace(-2,1.5,20); y = linspace(-2,2,20); [X,Y] = meshgrid(x,y); Now that we've created our grid of points, we can calculate the slope at each location using our previously-defined dydx function to dy represent == = ! p xX+y S 1l dydx(X,Y) m = 26x20 -4.0000 -3.8158 -3.6316 -3.4474 -3.2632 -3.0789 -2.8947 -2.7105 -2.5263 -2.3421 -2.1579 -1.9737 -1. -3.7895 -3.6053 -3.4211 -3.2368 -3.0526 -2.8684 -2.6842 -2.5000 -2.3158 -2.1316 -1.9474 -1.7632 -1. -3.5789 -3.3947 -3.2105 -3.0263 -2.8421 -2.6579 -2.4737 -2.2895 -2.1053 -1.9211 -1.7368 -1.5526 -1. -3.3684 -3.1842 -3.0000 -2.8158 -2.6316 -2.4474 -2.2632 -2.0789 -1.8947 -1.7105 -1.5263 -1.3421 -1. -3.1579 -2.9737 -2.7895 -2.6053 -2.4211 -2.2368 -2.0526 -1.8684 -1.6842 -1.5000 -1.3158 -1.1316 -0. -2.9474 -2.7632 -2.5789 -2.3947 -2.21605 -2.0263 -1.8421 -1.6579 -1.4737 -1.2895 -1.1053 -0.9211 -0. -2.7368 -2.5526 -2.3684 -2.1842 -2.0000 -1.8158 -1.6316 -1.4474 -1.2632 -1.0789 -0.8947 -0.7105 -0. -2.5263 -2.3421 -2.1579 -1.9737 -1.7895 -1.6053 -1.4211 -1.2368 -1.0526 -0.8684 -0.6842 -0.5000 -0. -2.3158 -2.1316 -1.9474 -1.7632 -1.5789 -1.3947 -1.21065 -1.0263 -0.8421 -0.6579 -0.4737 -0.2895 -0. -2.1053 -1.9211 -1.7368 -1.5526 -1.3684 -1.1842 -1.0000 -0.8158 -0.6316 -0.4474 -0.2632 -0.0789 Note that X, Y, and m are all now a 20 x 20 arrays. X and Y contain the coordinates of each (x, y) point, and m contains the of the slope of the function y(x) at each point. To use the quiver function to plot the slope at each point as an arrow, we need to break m up into its horizontal and vertical components, which we can do using the definition of slope: v m=— u Where 1 = Ax and v = Ay represent the "run" and "rise" of the slope at each (x, y) point. It follows then that:
Since quiver will scale the size of the arrows representing the slope m, we can assume that 4 = 1 everywhere, and then calculate v from m and wu. u = ones(size(X)) u = 20x20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 vV = m.*u vV = 20x20 -4.0000 -3.8158 -3.7895 -3.6053 -3.5789 -3.3947 -3.3684 -3.1842 -3.1579 -2.9737 -2.9474 -2.7632 -2.7368 -2.5526 -2.5263 -2.3421 -2.3158 -2.1316 -2.1053 -1.9211 Again, notice that u and v have the same dimensions as X, Y, and m. Also notice that the element-wise multiplication . * is used to ensure that each element in « is multiplied by the corresponding slope in m. Now u, v, X, andY are all defined, we are ready to plot using quiver. R R R R R RRRRR -3. 4211 -3. -3. -2. -2. -2. -2. -1. -1. -3 R R PR RPRRRRPRR 6316 2105 0000 7895 5789 3684 1579 9474 7368 R R R R R R RRBRRR@ 4474 .2368 .0263 .8158 .6053 .3947 .1842 .9737 .7632 .5526 quiver(X,Y,u,v, 'linewidth',1.25); grid on; ax = gca; PR R PR RPRRPRRRLR RR -3. -3. -2. -2 -2 R R PR RRRRRPRER 2632 0526 8421 .6316 4211 -2. -2. -1. -1. -1. 2105 0000 7895 5789 3684 R R PR RPRRRRPRPR vV =mu R R R R R RRRRR .0789 . 8684 .6579 4474 .2368 .0263 .8158 .6053 .3947 .1842 R R R R R RRRRR -2. -2. 4737 .2632 -2. -1. -1. 4211 -2 -2 -1 -1. -1. 8947 6842 0526 8421 6316 2105 0000 xlabel('x"'); ylabel('y'); title('Activity 3: Slope Field'); ax.FontSize = 12; xlim([-2,1]); ylim([-2.25,2.25]); ax.XAxisLocation = 'origin'; ax.YAxisLocation = ‘origin’; R R PR RPRRRRPRR R R R R R R RRBRRR@ . 7105 .5000 .2895 .0789 .8684 .6579 4474 .2368 .0263 .8158 P R PR RPRRRRLR PR R R R R RRRRRER .5263 .3158 .1053 .8947 .6842 4737 .2632 .0526 .8421 .6316 R R PR RPRRRRPRPR .3421 .1316 .9211 .7105 .5000 .2895 .0789 .8684 .6579 4474 R R R R R R RRRAR R R R R R RRRRR .1579 .9474 .7368 .5263 .3158 .1053 .8947 .6842 4737 .2632 R R PR RPRRRRPRR R R PR RPRRPRRRPRPR .9737 .7632 .5526 .3421 .1316 .9211 .7105 .5000 .2895 .0789
Activity 3: Slope Field \ I o I\ VA v A L A A I i | | i \ NN \\\01\\\ ¥ VY NN | | | \ | e 3 IR B I | \ ' ¥V NMNYNYD LAY M NN a | \ Y EIIEEEEREEET | | | (&) A28 e 8 1y vk N Y ANNYNNNNNN MAVNY Y N NNNYNYNN | LA AAAAAAT LAY NI NN NN NN I wn \ P4V 47 & A N di 4 i E A Al VA A | (&) S Xd Nt o nl g Iy s s agseoh | E Ik | A | | AT ANy AR AT S L4 £|4 ] A AR il o AN A A i q & 2 2 &l " Ly yyy £~V 7 dd It Ji P I P PPl PEEINN veveds L e s 200 108 eIV P W AT R N Ay e dd &ie 2E L 1 Pl LS AL TN P F AT e R A 4.7 Now that we’ve visualized what our slope field look like, and we have our ODE all set up, we can start using MATLAB's built-in functions to solve our ODEs. MATLAB's Non-Stiff ODE Solvers: ode23 and ode45 In Module 4, we introduced the Runge-Kutta techniques for numerically solving ODEs. Specifically, we showed that this family of techniques uses finite-difference approximations of the derivative at a point f'(x;.y;) to estimate the value of f(x;;1,yi+1). By increasing the order of the finite-difference approximation, we could increase the accuracy of our solution. While we now know how to write our own algorithms for performing these solution techniques, MATLAB already has a number of built-in Runge-Kutta solvers that we can use. For non-stiff ODEs, the most commonly used solvers are ode23, which performs the second-order Runge-Kutta method, and ode45, which performs the fourth-order Runge-Kutta method. Both functions use the same input argument structure. [t,y] = ode23(odefun,tspan,y0); [t,y] ode45(odefun,tspan,y0); Where the input variable tspan contains the range of independent variable values over which the ODE should be solved (usually a vector of x or y values), and y@ contains the initial condition y(0) = y,. The odefun input variable is a function handle referencing the ODE presented in standard form. In other words, whatever expression is contained within the odefun function must have the form: dy _ dx_f(x’y) The ODE function dydx in Activity 1 is in standard form. The output variables t and y are vectors containing the values of the solution y(#) and the independent variable ¢, respectively. Activity 4: Solving with ode23 and ode45 Let's see how we can now solve our ODE from Activity 3 using ode23 and ode45. Since we already have our differential equation set up and ready, we can simply use the solvers. % Select the initialpoint (x0,y0) X0 = -1.6; yo = 0.95; % Solve the ODE [x23,y23] = ode23(dydx,[x0,1],y0);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
[x45,y45] = ode45(dydx,[x0,1],y0); % Plot the solutions figure(l) hold on plot(x0,y0, 'o', 'linewidth',1.5); plot(x23,y23,'--', " "linewidth',1.5); plot(x45,y45,"'-.", "linewidth',1.5); hold off s = sprintf('(%1.2f,%1.2f)",[x0,y0]); title('Activity 4: ode23 and ode45 Solutions'); legend('Slope Field',s, ‘ode23','oded45', 'location’, 'eastoutside’); Activity 4: ode23 and ode45 Solutions I L4 P F a8 B /7 / :Q /,'#/ f PR e o b - @ W PR /',’/l/‘./'//,/,l \s-s-—-—--“;o,/'///"’:]//'/’/’/ - - ! K e S o"‘:-‘--u—: - /’4////"// )/ {/ /7 . 1. E:;Z‘:’fl_:—fl“' /’.«"'/'1/’/’ 4 7 ~ ~~~.______,fl' -, 2 2 N SRFIERE T 2 S W W . e A e o . N N R N SeRIsETE - =T 24 | T Slope Field s\. N - ’/ A £ % SN N T [ aTe e O (-1.60,0.95) NN BN Ny s s s == 7 7 e~ —ode23 B AAR NN N - ‘0’5‘ = * 1 |=-=-=0de45 ! TR - - MR SRR R VR T T R ; kALY M Y eSS e o B & A D R ' \ . w4 TH T e V30N N % N b N N B4R A EE\ TR ERR AR \\\\\‘\i‘q\\\\\\\\_{\\\\\ \\\\\\\\\\\\\\‘SZ\\\\w Use the slider-bars to change the initial point y(xo) = yo and see how the solution to the ODE changes. For this particular ODE, there is very little difference between the ode23 and ode45 solutions, but this will not always be the case. Activity 5: User-Define Step Size. In Activity 4, notice that changing the initial condition (xo, yo) sometimes changes the length of the results vectors x23, y23, x45, y45. This is because MATLAB will automatically select the step size for the Runge-Kutta algorithms based on the tspan input argument. If, however, you want to manually specify the Runge-Kutta step size, 4, you can do this by replacing the tspan range with a vector of points. To illustrate, let's solve the ODE: Z—x=—2x over the range —2 < x < 2 with initial condition y(—2) = —1. We'll first let MATLAB select the step size, and then use we'll solve with a user-selected step-size and compare the results: % Define the ODE to be solved odefun = @(x,y) -2*x; % Let MATLAB select h [ XAuto,yAuto] = ode23(odefun,[-2,2],-1)
XAuto = 13x1 -2.0000 -1.9800 -1.8800 -1.4800 -1.0800 -0.6800 -0.2800 0.1200 0.5200 0.9200 yAuto = 13x1 -1.0000 -0.9204 -0.5344 ©.8096 1.8336 2.5376 2.9216 2,9856 2.7296 2.1536 % User-Selected h h =#0.15 h = 9.1500 [xUser,yUser] = ode23(odefun,-2:h:2,-1) XUser = 27x1 -2.0000 -1.8500 -1.7000 -1.5500 -1.4000 -1.2500 -1.1000 -0.9500 -0.8000 -0.6500 yUser = 27x1 -1.0000 -90.4225 .1100 .5975 .0400 .4375 .7900 .0975 .3600 .5775 N NN P PP OO % Plot figure; plot(linspace(-2,2),3-1linspace(-2,2).72, "'k--', 'linewidth’,1.25); hold on plot(xAuto,yAuto, ‘ro-"', "'markerfacecolor', 'r'); plot(xUser,yUser, 'bd-', "markerfacecolor’,'b"); hold off grid on;
L = legend('Analytical Solution: $y = 3-x72%$', 'MATLAB-selected h', 'User-Selected h'); L.Interpreter = 'latex’'; L.Location = 'south’; xlabel('x"); ylabel('y'); 3 T T Y Wy 1 1 H 7 25} ) ' 5| S 15} | >~ 1} 0.5} | ol & = == = Analytical Solution: y =3 X 05T —&@— MATLAB-selected h ] —4@— User-Selected h -1 1 1 1 1 L : 1 -2 1.5 -1 0.5 0 0.5 1 1.5 2 X Solving Systems of ODEs MATLAB’s ODE solvers, like ode23 and ode45, can also work with function handles containing systems of equations. Recall from lecture, that a system of ODEs has the form: —fl(-xs _Yl, _)’2’ cee ,,Yn)- Y, (X, y1,y2, -2 s Yn) _fn(-x’ yl )’2, cee oy ,Yn)_ Where each row contains an expression for the ODE to be solved. Each of these expressions is a function of the independent variable , as well as one or more of the solution functions y;(x), y2(x), ... Note that each row represents an ODE in standard form. Solving this system of equations means finding the solution set: [y;(x), y2(x), ..., y.(x)] that satisfies all of the ODEs in the system. When using ode45, we can incorporate systems of equations by defining a function that contains a vector of expressions. We can illustrate with an example. Activity 6: Solving a System of ODEs Consider the following system of equations: dy = -=0.5 dx Y1 yi(0) = 4 dy> y»(0) = —= = 4-03y;—-0.1 , y2 Vi Which we want to solve over the range () < x < 2 with a step size of 1 = 0.5 To use ode45, we first need to represent this system of equations as an anonymous function. This time, rather than single equation, this function must now represent a vector of equations.
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
odeSystem = @(x,y) [-0.5*%y(1);4-0.3*y(2) - 0.1*y(1)]; Notice that the odeSystem function contains a 2 x | array of equations. The first row represents the equation for % and the second row represents the equation for % Similarly, the input argument y must now be a two-element array, where y(1) represents the current value of y,(x;), and y(2) represents the current value of y,(x;). Even though neither equation in this specific system depends directly on x, it is still required as an input argument. Since we have two equations, we also need to include two initial conditions, which we will define as a single column vector: yo = [4;6] yo = 2x1 4 6 Finally, we need to define a vector x to represent the range over which we want to solve the ODE: X = 0:0.5:2; Now we are ready to solve the ODE. Let's solve with ode23. [~,y23] = ode23(odeSystem,x,y0) y23 = 5x2 4.0000 6.0000 3.1152 6.8577 2.4261 7.6321 1.8894 8.3269 1.4715 8.9468 Notice that we do not need the x output argument since it will be identical to the input argument x. Looking at our results, we see that the y23 and y45 variables have been created, and they each contain two columns, which correspond to our solutions [y;y,], because we specified that y(1) = y; and y(2) = y3 in the odeSystem function. Let’s plot our results for y; and y, for both ode23 and ode45 . figure; plot(x,y23, 'o-"', 'linewidth’',1.5); xlabel('x"); ylabel('y'); ax = gca; title('Activity 6'); grid on; ylim([0,10]); L = legend('$y_1(x)$", "$y_2(x)$"); L.Location = "northwest’; L.Interpreter = "latex’'; L.FontSize = 11;
i ' Actl\'nty 6 Solving Second-Order ODEs Since they rely on numerical algorithms like the Runge-Kutta method, MATLAB’s ODE solvers can only operate on function handles containing first-order ODEs in standard form. Luckily, as we discussed in class, higher-order ODEs can be re-expressed as first-order ODEs using the Standard form. As an example, let’s consider the Van der Pol equation: L . dx d*x L _ Using x and x instead of = and e for simplicity, we have: t ¥—pu(l=x)i+x=0 Conversion to canonical form is based on the definition x = %()'c), and x = %(x). So letting y; = x, we can define the following equations: yoo= x 2 = x=y, y3 = X=y, Rearranging this set of equations into standard form produces: y, = » Y, = ¥3 Rewriting our original equation in terms of y,, y,, and y; produces: ys=p(1=y1)y2—y Which we can substitute into our system of standard-form ODEs to produce:
Y1 Y2 Yol (1= y7)y2 =y Now with this form, we are ready to solve with ode45! All we need is a range over which to solve, and some initial conditions. Activity 7 Let’s solve the Van der Pol equation for initial conditions: x(0) = 2 and x(0) = 0 over the range 0 <t < 20, with x =1 First, let’'s define our set of ODEs: mu = 1; vdpODE = @(t,y) [y(2); mu*(1-y(1)*2)*y(2) - y(1)] % ODE system representing Van der Pol equation vdpODE = function_handle with value: @(t,y) [y(2);mu*(1-y(1)"2)*y(2)-y(1)] Now we can define initial conditions and our range: yol = 1; ye2 = -1; yo = [y0l;y02]; % Initial conditions [x(@), xdot(®)] = [y1(©), y2(0)] t = linspace(0,20,100); And we are all set to solve the ODE: [tNew,y] = ode45(vdpODE,t,y0) tNew = 1606x1 %) .2020 .4040 .6061 .8081 .0101 .2121 4141 .6162 .8182 R R R RROOOO® 1.0000 -1.0000 0.7760 -1.2266 0.5001 -1.5183 0.1571 -1.8888 -0.2657 -2.2958 -0.7600 -2.5445 -1.2617 -2.3110 -1.6603 -1.5799 -1.8925 -0.7404 -1.9747 -0.1223 Once again, we see that y has two columns, corresponding to [y;, y»],. Remembering that we defined y, = x and y, = x, we know that these columns really correspond to our solution x and its derivative: [y, y2] = [x, x] % First column of y refers to solution x X = y(:Jl); % Second column of y refers to the derivative dx/dt xdot = y(:,2); % Plot plot(tNew, X, tNew,xdot, 'linewidth',1.5); xlabel('t (s)'); ylabel('x'); ax = gca; title('Van der Pol Results'); ax.FontSize = 11; grid on; L = legend('x(t)", "$\dot{x}(t)$"', "location’, "'northwest’'); L.FontSize = 11; L.Interpreter = "latex’;
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
Van der Pol Results x(t) x(1) -3 1 1 L 0 5 10 15 20 t(s) Post Lab Complete the Lab 12 Post Lab Assignment on MATLAB Grader References and Documentation [1] https://www.mathworks.com/help/matlab/ref/meshgrid.html [2] https://www.mathworks.com/help/matlab/ref/ode23.html [3] https://www.mathworks.com/help/matlab/ref/ode45.html