20_MATLAB_ode45W12

20_MATLAB_ode45W12 - Page - 1 - of 10 Numerical Method to...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Page - 1 - of 10 Numerical Method to Solve Initial Value Problems (IVPs) CHE 361 Both Euler and RK4 use fixed step size in time. (Can be used with nonlinear as well as linear ordinary differential Equations, ODEs) __________________________________________________________________________ Given either a single equation or a set of first-order differential equations of the form: dxi (t ) = f= RHSi , with i = 1, ... , n (for a set of ODEs) i [ x (t ), u (t ) ] dt Start = t0 , with x x0= u= u0 , calculate each xi at one "step" (∆t ) later in time. at t = ,u (t0 ) __________________________________________________________________________ 1. Euler Method (similar to Eq 7-29, pg 125 SEMD3) Simplest method, very small steps usually required. Concept is that “average slope” of x vs. t is the “initial” slope. xi (t0 + ∆t )= xi (t0 ) + { fi [ x(t0 ) , u (t0 ) ]}∆t= xi + RHSi (t0 ) ∆t 2. Runga-Kutta Method (classical 4th-order RK4 = see also Sec. 2.5 pg 33 SEMD3 discussion) xi (t0 + ∆t ) =t0 ) + xi ( 1 xi {Ki,A + 2 Ki,B + 2 Ki,C + Ki,D }∆t =+ Ki, average ∆t 6 dx Each Ki,J can be considered to be an estimate of the average rate of change dt during the time interval ∆t . The calculated “average slope” ( Ki, average ) is a weighted sum of the 4 “estimates”: K i ,A = f= i [ x 0 , u (t0 ) ] RHSi (t0 ) K A ∆t ∆t K i ,B = fi x 0 + , u (t0 + ) 2 2 K B ∆t ∆t K i ,C = fi x 0 + , u (t0 + ) 2 2 K i ,D = fi x 0 + K C ∆t , u (t0 + ∆t ) 3. Runga-Kutta(Dormand-Prince) = MATLAB ode45 and other explicit variable ∆t time step methods. These are explicit in the sense that calculations do not involve solutions of implicit equations at each time step. The method automatically adjusts the time step to enforce specified tolerances on the estimated error. “45” = 4th-order accuracy using a 5th-order check (to get same result at next time point). Page - 2 - of 10 SIMPLE EXAMPLE : Euler Method for single ODE dy (t ) = t ) + u (t ) + 2 − y( dt with u =4 for t < 0 and u (t ) = t ≥ 0. 2, y = 2.2 for At t = 0, y = 4, dy/dt = -4 + 2.2 + 2 = 0.2 Now for a time step of Δt = 0.5, y(0.5) = y(0) + (dy/dt at 0)*Δt = 4 + 0.2*0.5 = 4.1 Then (dy/dt at 0.5) = -y(0.5) + 4.2 = -4.1 + 4.2 = 0.1 and y(1) = y(0.5) + (dy/dt at 0.5)*Δt = 4.1 + 0.1*0.5 = 4.15 Then (dy/dt at 1) = -y(1) + 4.2 = -4.15 + 4.2 = 0.05 and y(1.5) = y(1) + (dy/dt at 1)*Δt = 4.15 + 0.05*0.5 = 4.175 ...and so forth until dy/dt = 0 at the new steady state Using Laplace transforms, we can easily find the analytical (exact) solution to this linear ODE and evaluate the accuracy of the Euler method: y (t )analytical 4.2 − 0.2e − t = Time yanalytical = Exact yEuler, Δt = 0.5 dy/dt = 4.2 - yEuler 0 4.0000 4.000 0.20 0.5 4.0787 4.100 0.10 1. 4.1264 4.150 0.05 1.5 4.1554 4.175 0.025 2. 4.1729 4.1875 0.0125 2.5 4.1836 4.19375 0.00625 Improved accuracy over a time period can be obtained by taking more, but smaller steps. For extremely small steps, the larger number of steps may be slow and roundoff errors may lead to less accurate results. Page - 3 - of 10 COMPLEX EXAMPLE : Euler and RK4 methods for 2 ODEs for 2 states: x1 and x2 dx1 (t ) = x2 (t ) dt dx2 (t ) = x1 (t ) − 3 x2 (t ) + 12.5u (t ) −6.25 dt with [u = 1 = 2 =for t < 0 and u = t ≥ 0. 0, x 0, x 0] 1 for It can be shown the the solution for x1(t) would be the same as for y(t) for a unit step to the transfer function 2 Y (s) G ( s) , 2nd-order with τ =0.4, ζ =0.6 = = 2 0.16 s + 0.48s + 1 U ( s ) The analytical solution using Laplace transforms then can be used for comparisons: 2{ 1 − e y (t )analytical =−1.5t [ cos(2t ) + 0.75sin(2t ) ] } For example at t = 1 we get the exact answer to be yanalytical = 1.881372, compared to yEuler = 1.964 with a time step of 0.05. For a time step 4 times longer (0.2 vs. 0.05), the RK4 method yields yRK4 = 1.882, which is more accurate than Euler’s result. The increased complexity of RK4 allows larger steps in time with better accuracy. The following Tables were done using EXCE spreadsheet software. The Simulink block diagram num-meth.mdl was used to draw plots of numerical methods vs. exact solutions and show similar results as found in the tables. t 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 x1=y 0.000 0.000 0.031 0.089 0.169 0.267 0.378 0.501 0.630 0.763 0.897 1.031 1.162 1.288 1.409 1.522 1.628 1.726 1.814 1.894 1.964 x2 0.000 0.625 1.156 1.598 1.956 2.234 2.441 2.581 2.663 2.692 2.675 2.618 2.528 2.411 2.272 2.116 1.947 1.772 1.592 1.411 1.232 dx1/dt 0.000 0.625 1.156 1.598 1.956 2.234 2.441 2.581 2.663 2.692 2.675 2.618 2.528 2.411 2.272 2.116 1.947 1.772 1.592 1.411 1.232 dx2/dt 12.500 10.625 8.836 7.149 5.577 4.130 2.812 1.627 0.577 -0.342 -1.132 -1.798 -2.346 -2.784 -3.120 -3.362 -3.519 -3.600 -3.613 -3.569 -3.474 delt 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 Page - 4 - of 10 Results from RK4 method via Excel for complex example with a time step of Δt = 0.2 t 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 x1=y 0.000 0.2023 0.645 1.137 1.567 1.882 2.077 2.170 2.190 2.165 2.122 2.075 2.036 2.008 1.991 1.983 1.982 1.985 1.989 1.993 1.997 x2 0.000 1.805 2.463 2.370 1.882 1.268 0.697 0.255 -0.035 -0.187 -0.236 -0.220 -0.170 -0.112 -0.059 -0.019 0.006 0.019 0.023 0.020 0.015 dx1/dt 0.000 1.805 2.463 2.370 1.882 1.268 0.697 0.255 -0.035 -0.187 -0.236 -0.220 -0.170 -0.112 -0.059 -0.019 0.006 0.019 0.023 0.020 0.015 dx2/dt 12.500 5.822 1.083 -1.717 -2.938 -3.067 -2.572 -1.827 -1.081 -0.472 -0.051 0.189 0.287 0.287 0.234 0.162 0.093 0.038 0.000 -0.020 -0.028 delt 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 K1,A 0.000 1.805 2.463 2.370 1.882 1.268 0.697 0.255 -0.035 -0.187 -0.236 -0.220 -0.170 -0.112 -0.059 -0.019 0.006 0.019 0.023 0.020 0.015 K2,A 12.500 5.822 1.083 -1.717 -2.938 -3.067 -2.572 -1.827 -1.081 -0.472 -0.051 0.189 0.287 0.287 0.234 0.162 0.093 0.038 0.000 -0.020 -0.028 K1,B 1.250 2.387 2.571 2.198 1.588 0.961 0.439 0.072 -0.143 -0.234 -0.241 -0.201 -0.142 -0.083 -0.035 -0.003 0.016 0.023 0.023 0.018 0.013 K2,B 8.750 2.947 -0.781 -2.683 -3.233 -2.939 -2.236 -1.438 -0.735 -0.214 0.112 0.270 0.307 0.271 0.201 0.126 0.061 0.014 -0.014 -0.027 -0.029 K1,C 0.875 2.099 2.384 2.102 1.559 0.974 0.473 0.111 -0.108 -0.208 -0.225 -0.193 -0.140 -0.085 -0.039 -0.007 0.012 0.020 0.021 0.018 0.012 K2,C 9.094 3.446 -0.289 -2.286 -2.961 -2.786 -2.176 -1.440 -0.771 -0.262 0.066 0.234 0.283 0.258 0.196 0.126 0.065 0.019 -0.010 -0.024 -0.027 K1,D 1.819 2.494 2.405 1.913 1.290 0.711 0.261 -0.033 -0.189 -0.239 -0.223 -0.173 -0.114 -0.060 -0.020 0.006 0.019 0.023 0.021 0.016 0.010 K2,D 5.950 1.130 -1.724 -2.973 -3.110 -2.613 -1.858 -1.101 -0.483 -0.055 0.191 0.290 0.291 0.238 0.165 0.095 0.038 0.001 -0.020 -0.028 -0.027 K1avg 1.0115 2.212 2.463 2.147 1.578 0.975 0.464 0.098 -0.121 -0.219 -0.232 -0.197 -0.141 -0.084 -0.038 -0.005 0.014 0.021 0.022 0.018 0.013 K2avg 9.023 3.290 -0.463 -2.438 -3.073 -2.855 -2.209 -1.447 -0.763 -0.246 0.083 0.248 0.293 0.264 0.199 0.127 0.064 0.018 -0.011 -0.025 -0.028 Page - 5 - of 10 num_meth.mdl "complex" example = 2 odes 2 Step_Inout 0.16s2 +0.48s+1 TF_model Mux ErrorPlot Clock Interpreted MATLAB Fcn y_exact Numerical Method (line) vs. Exact Solution (circles) for Constant Time Steps 2.5 Output, y-NumMeth = line, y-analytical = o Mux 2 Euler Method with Delta t = 0.2 1.5 1 0.5 0 0 0.2 0.4 0.6 0.8 1 Time , t 2 Euler Method with Delta t = 0.05 1.5 1 0.5 0 0 0.2 0.4 0.6 Time , t Numerical Method (line) vs. Exact Solution (circles) for Constant Time Steps 2.5 Output, y-NumMeth = line, y-analytical = o Output, y-NumMeth = line, y-analytical = o Numerical Method (line) vs. Exact Solution (circles) for Constant Time Steps 2.5 0.8 1 2 RK4 Method with Delta t = 0.2 1.5 1 0.5 0 0 0.2 0.6 0.4 Time , t 0.8 1 Page - 6 - of 10 Solving Linear ODEs using MATLAB = ode45 function (CHE 361) Here we solve a linear ODEs example of using the MATLAB built-in function ode45. We need to write 2 files: Ex2_odes.m and RHS2.m below. ode45 also works with nonlinear ODEs. dx1 (t ) = x2 (t ) and dt dx2 (t ) = x1 (t ) − 3 x2 (t ) + 12.5u (t ) with u (t ) = t ≥ 0 with x1 = at t 0. −6.25 1 for x2 == 0 dt % Ex2_odes.m – Dr. K. Levien OSU – CHE 361 clear all, format compact yinitial = [0 0]; tstop = input('Enter simulation stop time in hours (1) : '); tspan = [0:0.05:tstop]; % calculate values at increments of 0.05 [t,y] = ode45(@RHS2,tspan,yinitial); x1 = y(:,1); % x1 is first column in y matrix x2 = y(:,2); % x2 is second column in y matrix figure plot(t,x1,'-b',t,x2,'r--') % x1 as solid blue line, x2 as dashed red line. grid % puts a grid on plot to help read values xlabel('Time') ylabel('x_1 = solid, x_2 = dashed') title('Example with 2 ODEs from Numerical Methods Handout: CHE 361') Results = [t';y']' % Prints out Results as t,x1,x2 in columns function dydt = RHS2(t,y) x1 = y(1,:); x2 = y(2,:); dydt(1,1) = x2; dydt(2,1) = -6.25*x1 - 3*x2 + 12.5; % Here u(t) = 1 for 12.5u(t) >> Ex2_odes Enter simulation stop time in hours (1) : 1 Results = Example with 2 ODEs from Numerical Methods Handout: CHE 361 0 0.0149 0.0564 0.1203 0.2026 0.2994 0.4074 0.5235 0.6447 0.7688 0.8933 1.0166 1.1369 1.2530 1.3638 1.4683 1.5660 1.6564 1.7391 1.8141 1.8814 0 0.5789 1.0687 1.4749 1.8030 2.0594 2.2502 2.3818 2.4606 2.4927 2.4843 2.4410 2.3684 2.2715 2.1553 2.0240 1.8817 1.7319 1.5779 1.4225 1.2681 2.5 2 x1 = solid, x2 = dashed 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.5 1 0.5 0 -0.5 0 0.5 1 1.5 2 2.5 3 Time NOTE: The final answer using ode45 here is x1 = 1.8814, which is the same to 5 significant digits as the analytical result in handout of x1 = 1.881372 at time = 1.0. Plot made for Time up to 3 to show overshoot. The Figure Editor was used to increase font size/line widths and use a solid line grid. Page - 7 - of 10 Solving Nonlinear ODEs using MATLAB Use of “ode45” for bioreactor ODEs (CHE 361 - see Appendix C, SEM2/SEMD3) The main “driver” program = bioreactor45.m % bioreactor45.m clear all, format compact % Simulation of the CHE 361 Bioreactor project as an initial value problem % for 2 nonlinear ODEs. Dr. K. Levien, Oregon State Univ., Chemical Engineering % help ode45 %[T,Y] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS,P1,P2...) passes the additional % parameters P1,P2,... to the ODE function as ODEFUN(T,Y,P1,P2...), and to % all functions specified in OPTIONS. Use OPTIONS = as a place holder if % no options are set. % First define values for the 4 kinetic parameters and volume of the reactor. % These values are for the base case bioreactor. M = 0.1; mumax = 2.; Y = 0.5; KS = 1.; V = 1; % One liter = reactor volume % Next define the inputs, which are constant for t greater than or equal to 0. F = 1.1; % Here flow is 1.1 liters/hr. The nominal steady state has F = 1.0 Ni = 5.2; % Nominal steady-state nutrient feed concentration is 5.2 gms/liter. % The initial values here are those corresponding to the nominal steady state, % but they could have any reasonable and realistic values. If Binitial = 0, % for example, there will never be any bugs in the reactor, B = 0 always. Binitial = 2; Ninitial = 1; % This m-file is written to be interactive: it asks for a value of the % final time to be entered in the MATLAB Command Window. tstop = input('Enter simulation stop time in hours (8) : '); tspan = [0 tstop]; % % % % % % % % % % The list of arguments to ode45 starts with the name of the function. Here the user supplied function must be called bioRHSs.m and follow ] the rules/syntax for a MATLAB function. The next vector contains at least an initial value of time and the desired final value of time. The next vector contains the initial values of the “states”. The next argument is the null vector, which indicates the default options for the ode45 function are to be used. Following the placeholder for the “options” vector are the individual constants which are needed in the RHS equations. These can be scalers or vectors. Here they are all scalers. [t,y] = ode45(@bioRHSs,tspan,[Binitial,Ninitial],,M,mumax,Y,KS,V,F,Ni); % The output from ode45 is a vector of time values and a matrix of y values. % The columns of the y matrix are the vectors of values of the outputs. % Here the first output (column) is B, the second output (column) is N. B = y(:,1); N = y(:,2); Page - 8 - of 10 % Now that we have calculated the solution of the two ODEs, we plot each % output as a function of time in a separate Figure Window with appropriate % labels for the x and y axes and t title to appear at the top of the plot. figure plot(t,B,'-bo','LineWidth',1.5) xlabel('Time , h') ylabel('Product Bug Concentration , g/L') title('Dynamic Simulation of Bioreactor, B(t) from ode45 in MATLAB') figure plot(t,N,'--r*','LineWidth',1.5) xlabel('Time , h') ylabel('Product Nutrient Concentration , g/L') title('Dynamic Simulation of Bioreactor, N(t) from ode45 in MATLAB') % The “user-supplied” function file to calculate the two “right hand sides” % (RHSs) of the bioreactor first-order ordinary differential equations (ODEs) function dydt = bioRHSs(t,y,M,mumax,Y,KS,V,F,Ni) B = y(1); % The 1st “output” in the y output vector is B(t). dy(1,1)/dt = dB/dt N = y(2); % The 2nd “output” in the y output vector is N(t). dy(2,1)/dt = dN/dt mu = mumax*(N/(KS + N)); % This calculates the “intermediate” variable mu. N KS + N µ = µmax % The dB/dt expression, Eqn (2) of the handout. dB FB F = µB − = B µ − dt V V dydt(1,1) = B*(mu – F/V); % The dN/dt expression, Eqn (4) of the handout. ( µ B ) − MB dN FN i FN µ B F = − − − MB = ( N i − N ) − dt V V Y Y V dydt(2,1) = (F/V)*(Ni - N) - (mu*B)/Y – M*B; ------------------------------------------------------- Page - 9 - of 10 The interactions in the MATLAB Command Window look like: >> bioreactor45 Enter simulation stop time in hours (8) : 8. The following command displays results in three columns: t, B, N t is a column vector, y is a matrix with 2 collumns. 1.0000 1.0447 1.0797 1.1075 1.1297 1.1519 1.1689 1.1818 1.1917 1.2014 1.2083 1.2130 1.2162 1.2186 1.2202 1.2212 1.2218 1.2223 1.2225 1.2226 1.2227 1.2227 1.2227 1.2226 1.2226 1.2225 1.2225 1.2225 1.2224 1.2224 1.2224 1.2223 1.2223 1.2223 1.2223 1.2223 1.2223 1.2223 1.2222 1.2222 1.2222 1.2222 1.2222 1.2222 1.2222 Dynamic Simulation of Bioreactor, B(t) from ode45 in MATLAB 2 1.99 Product Bug Concentration , g/L 2.0000 1.9789 1.9625 1.9497 1.9396 1.9296 1.9222 1.9167 1.9126 1.9086 1.9059 1.9043 1.9032 1.9024 1.9020 1.9018 1.9017 1.9017 1.9017 1.9017 1.9018 1.9019 1.9020 1.9020 1.9021 1.9021 1.9022 1.9022 1.9023 1.9023 1.9023 1.9023 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.9024 1.98 1.97 1.96 1.95 1.94 1.93 1.92 1.91 1.9 0 1 2 3 4 Time , h 5 6 7 8 Dynamic Simulation of Bioreactor, N(t) from ode45 in MATLAB 1.25 Product Nutrient Concentration , g/L >> [t';y']' ans = 0 0.1196 0.2392 0.3588 0.4785 0.6288 0.7790 0.9293 1.0796 1.2762 1.4728 1.6693 1.8659 2.0659 2.2659 2.4659 2.6659 2.8659 3.0659 3.2659 3.4659 3.6659 3.8659 4.0659 4.2659 4.4659 4.6659 4.8659 5.0659 5.2659 5.4659 5.6659 5.8659 6.0659 6.2659 6.4659 6.6659 6.8659 7.0659 7.2659 7.4659 7.5994 7.7329 7.8665 8.0000 1.2 1.15 1.1 1.05 1 0 1 2 3 4 Time , h 5 6 7 8 Page - 10 - of 10 fnicefig.m = a utility function for plot modification Below is an example of a function that can be used to modify Simulink plots to have thicker lines (2), larger font size (12), and a font of Times New Roman, which is a common font for professional journals and publications. In addition grid lines are changed to solid lines. Note that this is a function file, not a program, so the standard “clear all” and “format compact” commands are not used. Workspace variables are thus retained so that intermediate or final values of variables can be examined. Since there is no output to the Command Window, the format compact command is not needed. The first command and last command for “ShowHiddenHandles” are needed when modifying plots made by Simulink using the standard CHE 361/461 plot block in a Simulink diagram. function fnicefig % Dr. K. Levien, Oregon State University % CHE 361 W2012 % A function to make a Simulink Figure "nicer". % Set size of axis tic labels and axes thickness. set(0,'ShowHiddenHandles','on') set(gca,'LineWidth',2,'FontName','Times New Roman','FontSize',12, ... 'GridLineStyle','-'), grid % Set all line thicknesses to 2. h1 = findobj(gca,'Type','line'); set(h1,'LineWidth',2) % Set xlabel, ylabel and h2=get(gca,'Xlabel'); set(h2,'FontName','Times % h3=get(gca,'Ylabel'); set(h3,'FontName','Times % h4=get(gca,'Title'); set(h4,'FontName','Times title font name and size. New Roman','FontSize',12) New Roman','FontSize',12) New Roman','FontSize',12) set(0,'ShowHiddenHandles','off') ...
View Full Document

Ask a homework question - tutors are online