Math 216: Differential EquationsLab 2: Eulerâs Method and RC CircuitsGoalsIn this lab you will implement Eulerâs method to approximate measurements of the chargeon a capacitor in a basic RC circuit. You will learn how to write .m files for Matlab andhow to program Eulerâs method; then you will investigate some of the limitations of Eulerâsmethod.Application: a basic RC circuitThe state of an electrical circuit consisting of resistors, capacitors, and an applied voltagecan be described by differential equations.Consider the following circuitISR+E(t)C-with resistance R Ohms (Ω), capacitance C Farads (F), and applied voltage E(t) Volts (V).The charge on the top plate of the capacitor at time t is Q(t) Coulombs (C), and the currentthrough the resistor is I(t) Amperes (A). The resistor has a voltage drop of RI, and thecapacitor has a voltage drop of Q/C. When switch S is closed at time t = 0, the sum of thevoltage across the resistor and the capacitor must equal the applied voltage. This gives usthe equation1RI + Q = E(t) .CThe current in the circuit is the rate of change of the amount of charge on the capacitor. Sousing the relationship dQ/dt = I, this becomes a first order differential equation for Q(t):RdQ1+ Q = E(t) .dtC1The initial condition for this equation is Q0 = Q(0), the initial amount of charge on thecapacitor. (Q0 could be set by imposing a voltage V0 = Q0 /C across the capacitor beforeinserting it into the circuit.) In this lab, we will use Eulerâs method to numerically solve thisdifferential equation for two different applied voltages: a constant voltage E0 , and then anAC voltage E(t) = 117 sin(120πt) which corresponds to the voltage out of a standard wallsocket in the United States.Prelab assignmentBefore arriving in lab, answer the following questions. You will need your answers in lab towork the problems, and your recitation instructor may check that you have brought them.These problems are to be handed in as part of your lab report.1. (a) We first consider the case where E(t) = E0 , a constant applied voltage. Verifythat the functionQ1 (t) = E0 C 1 − e−t/(RC)(1)is a solution to the initial value problem with E(t) = E0 ,11dQ1=−Q1 + E0 ,dtRCRQ1 (0) = 0 ,(2)where R, C, and E0 are constants. That is, by plugging t = 0 into the formula (1)show that the initial condition is satisfied, and then by differentiating the formula(1) and comparing with the right-hand side of the differential equation show thatQ1 (t) satisfies the differential equation. (In other words, do not try to find thesolution of the initial-value problem, but rather just check that the given functionsolves the problem.) Then use the exact solution formula (1) with R = 18000Ω,C = 0.0000125F,and E0 = 117V to complete the column labelled âExact yâ onTable 2 on the last page of the lab, for use in Lab problem 1.(b) Next consider the case where E(t) = E0 sin(120πt). Verify that the functionQ2 (t) =E0120π e−γt − cos(120πt) + γ sin(120πt)2 + (120π)2 )R (γ(3)satisfies the initial-value problem with this E(t),dQ211=−Q2 + E0 sin(120πt) ,dtRCRQ2 (0) = 0 ,(4)where the constant γ = 1/(RC) has units of Hz=sec−1 . What do the initial1condition and forcing term R E0 sin(120πt) represent for the system at the timethe switch is closed?22. Suppose you implement Eulerâs method using Matlab, using step size h, and createa vector t of time steps from t = 0 to t = 1. Often we refer to the first entry ast0 = 0, the next as t1 and the final entry will be tN = 1 where N h = 1. Matlab doesnot enumerate these entries in the same way. The first element of a vector in Matlabis always t(1). In this case, we will have t(1)=0, and t(N+1)=1. Find the Matlabindices n so that t(n)=0, t(n)=.5, t(n)=.86, and t(n)=1 if you used(a) N = 12 (Note: you cannot get t(n)=.86 in this case.)(b) N = 120 (Note: you cannot get t(n)=.86 in this case.)(c) N = 1200(d) N = 2400Record these values of n in Table 1 below.N1212012002400n so that t(n)=0n so that t(n)=.5n so that t(n)=.86n so that t(n)=1Table 1: Matlab Indices for Time VectorIn the labThe primary manner in which we will work with Matlab is by writing simple programs thatwe can then run in Matlab. Each program is a small file with a .m extension. For this lab wewill write programs that implement Eulerâs method, and will compare the results obtainedfor several step sizes h with exact solutions (when these are available).Creating a (new) fileLaunch Matlab. Use the buttons in the âCurrent Folderâ box to the left of the Matlabwindow to make sure that you are working on the Desktop. Note that we will work with fileson the Desktop in this labâbut that you will need to be sure to save them to your (UMich)IFS space, Google docs or Box.net (or Dropbox, etc.) account to have them for future useafter you finish in the lab today. Anything you do not move to one of these will be lost.Start editing a file by going to the File menu and selecting New and Script. A new windowcalled untitled will open; give it a name by selecting File→Save As in this window; give yourfile the name EULER.m and save it. You have now created an (empty) file called EULER.m.Note that capitalization matters here.3The contents of the EULER.m fileThis program will implement Eulerâs method to solve the initial-value problemdy= f (t, y) ,dty(a) = y0 .(5)In this application, the function y(t) stands for Q(t), the charge on the capacitor. In theprogram below, everything following a % is a comment. Comments give you informationabout the program, but are not evaluated by the computer. You may choose whether or notto type these comments into your program, but if you include the comments in your file youmust include the %, or the computer will try to read them. Type the following program intothe file EULER.m you have just created:clear t% Clears old time steps andclear y% y values from previous runsa=0;% Initial timeb=1;% Final timeN=12;% Number of time stepsy0=0;% Initial value y(a)h=(b-a)/N;% Time stept(1)=a;y(1)=y0;for n=1:N% For loop, sets next t,y valuest(n+1)=t(n)+h;y(n+1)=y(n)+h*f(t(n),y(n)); % Calls the function f(t,y)=dy/dtendplot(t,y)title([âEuler Method using N=â,num2str(N),â steps, by MYNAMEâ])% Include your own nameSubstitute your own name in place of MYNAME so it will appear on the plots you create. Whenyou are finished typing in the program statements, save your work by choosing Save fromthe File menu, or by clicking on the floppy disk icon on the EULER.m editor window.The right-hand side of the differential equation and the file f.mSince the program EULER.m refers to the function f(t,y), we must create a second .m filecalled f.m to define the right-hand side of the differential equation for Matlab. Go to the Filemenu and select New then Script. A new window called untitled will open. Choose Save As. . .from the File menu and a dialog box will open with a default entry untitled.m. Replacethat with f.m and click Save as you did before. You should then type into the file f.m thefollowing commands:4function f=f(t,y)R=18000;% ResistanceC=.0000125;% CapacitanceE0=117;% Constant Voltagef=-y/(R*C)+E0/R;% Defines the function fBe sure to save your work when you have finished typing.The file f.m contains the function f (t, y) for the right-hand side of the general differentialequation of the initial-value problem (5) above; the particular form of f (t, y) appearing on thelast line corresponds to the specific initial-value problem (2). To solve a different differentialequation with EULER.m or another solver, you need only change this file. The other lines inf.m serve to define specific values for the constants; here we are using R = 18kΩ, C = 12.5µF,and E0 = 117V.Note that your newly created f.m file gives the right-hand side of equation (5) for thecase of a constant forcing E(t): equation (2). We will change this later to solve in the caseof a non-constant forcing.The exact solution of the initial-value problem and the file yE.mSince for this problem we also know a formula for the exact solution to the initial-valueproblem, we can write a short Matlab program to evaluate this formula so we can comparewith the results of using Eulerâs method for various step sizes h. To tell Matlab about theexact solution formula, you should create another new file called yE.m. The file yE.m willcontain the exact solution Q1 (t) of the initial-value problem (2) corresponding to the righthand side function f (t, y) defined in the file f.m. If you solve a different differential equationwith EULER.m or one of the other numerical methods we will study in Lab 3, and you wishto compare with an analytical expression for the exact solution, you should modify the fileyE.m as well as f.m. Type the following commands into a new file called yE.m, saving yourwork when you have finished:function yE=yE(t)R=18000;% ResistanceC=.0000125;% CapacitanceE0=117;% Constant Voltagegam=1/(R*C);yE=E0*C*(1-exp(-gam*t)); % Exact solution yERunning your programSo far you have written a program and saved it. Now you want to use the program. In theCommand Window type EULER. The plot of the solution curve should then appear in a newwindow.5Here is a summary of the programming we have just done. To solve an initial-valueproblem for a differential equation you need a solver. In this case, EULER is your solver,and generates a numerical approximation to the solution to the differential equation usingEulerâs method. Then you must tell your solver what it is supposed to solve; that is yourf.m file defining the right-hand side of the differential equation. Finally, the solver must betold where to start, so you must specify initial conditions. In this case the initial conditionsare entered directly into the solver with the lines a=0 and y0=0. You âfine tuneâ your solverby changing N. The stopping time is also entered in the solver as b=1.Lab problems1. Solve initial-value problem (2) using EULER.m with N = 12.(a) Print your results as follows: After the plot has appeared, go to the File menu onthe figure window and choose Print. . . ; a dialog box will appear. If you click OKin the dialog box the plot will print on the printer in your lab.By the way, you can use the Edit menu on the figure window to change the sizeof your plot, add text, labels, legends, titles, etc.(b) Using your indices from Prelab problem 2 and your approximate solution vectory calculated by EULER.m, complete the portion of Table 2 corresponding to N =12. To view the nth component of the approximate solution vector y, just typey(n) in the Matlab command window and press return. Similarly, to view thecorresponding exact solution value, type yE(t(n)) and press return. Make sureyou have enough significant digits to compare the approximate and exact answers.You can see more significant digits in Matlab by typing format long in theMatlab command window and then hitting return.2. Repeat Problem 1 again using N = 120 and N = 1200. Use these data to completethe rest of Table 2.3. What guesses can you make about the error at a given time as N increases?4. Graphically investigate the error in solving the initial-value problem (4)âthe case ofa non-constant forcing function. To solve this problem using EULER.m, you will haveto change your file f.m so that the function value is different. Change the linef=-y/(R*C)+E0/R;to the linef=-y/(R*C)+(E0/R)*sin(120*pi*t);Note: probably the best way to do this is to comment the previous definition, but toleave it in the file, so that your file f.m becomes6function f=f(t,y)R=18000;C=.0000125;E0=117;% f=-y/(R*C)+E0/R;f=-y/(R*C)+(E0/R)*sin(120*pi*t);%%%%%ResistanceCapacitanceConstant VoltageOld definition of fNew definition of fTo compute the exact solution Q2 (t) for the initial-value problem (4), you will similarlyneed to modify the file yE.m. Change the lineyE=E0*C*(1-exp(-gam*t);to% yE=E0*C*(1-exp(-gam*t);And then add after that the linesomg=120*pi;A=E0*omg/(R*(gam^2+omg^2));B=E0*gam/(R*(gam^2+omg^2));yE=A*(exp(-gam*t)-cos(omg*t))+B*sin(omg*t);You can move among the various files in the editor window by clicking on the tabs atthe bottom of the window.To plot the approximate solution and the exact one on the same set of axes, replacethe line plot(t,y) near the end of the file EULER.m with the linesplot(t,y,â:â, t,yE(t),â-â)legend(âApproximateâ,âExactâ)This defines a vector corresponding to the exact solution. Sample values of the exactsolution are plotted connected by a solid line, and the Eulerâs method approximationis plotted on the same axes as a dotted line. The legend command puts a helpfulbox in the corner of the plot to help you identify which graph corresponds to whichfunction. You can change the location of the legend on the graph by clicking on it orby using the Edit menu on the figure.5. Run EULER using N = 120, N = 1200, and N = 2400. Be sure that you understandwhat is being graphed: for N = 120, we are plotting the Euler approximation to thesolution at the points t = 0, t = 0.008333, t = 0.016667, etc.âand are plotting thevalue of the exact solution, yE , at the same points. Explain how this results in thegraph of the exact solution being different when you use the different values of N .6. Many different things can give rise to error in numerical approximations. One is theaccuracy of the numerical method, which we saw when we considered the constantforcing functionâas we increased the number of steps, the accuracy of the approximation improved. Another, undersampling, comes in to play when an insufficient numberof points are used to graph the solution. Undersampling occurs when the data pointsare spaced too far apart to capture all the behavior of the equation. You may know7this phenomenon as âaliasing,â when an insufficiently sampled high frequency signalappears to be a lower frequency signal. Even if you have calculated the exact solutioncorrectly its graph may not be accurate. Imagine plotting a sine wave with frequency120 cycles per second by sampling the function 120 times a second; the sample pointswill suggest a constant function instead of an oscillatory function. Did you see thisproblem here? Try N = 50, N = 75 and N = 100. Do these result in aliasing?Note: your lab report for Lab 2 should consist of your solutions to the prelabproblems, your solutions to lab problems 1â6, and a brief, original, conclusionsparagraph summarizing what you have learned.TimeIndexExact yApproximate yErrortnyE(t)y(n)|y(n)−yE(t)|Percent Errory(n) − yE(t)à 100yE(t)N = 12.51N = 120.51N = 1200.5.861Table 2: Approximate Solutions to Equation (2)last modified 9 Sep 20158