Math 216: Differential Equations Lab 2: Euler’s Method and RC Circuits Goals

| August 31, 2017

Question
Math 216: Differential Equations
Lab 2: Euler’s Method and RC Circuits
Goals
In this lab you will implement Euler’s method to approximate measurements of the charge
on a capacitor in a basic RC circuit. You will learn how to write .m files for Matlab and
how to program Euler’s method; then you will investigate some of the limitations of Euler’s
method.

Application: a basic RC circuit
The state of an electrical circuit consisting of resistors, capacitors, and an applied voltage
can be described by differential equations.
Consider the following circuit

I
S

R

+

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 current
through the resistor is I(t) Amperes (A). The resistor has a voltage drop of RI, and the
capacitor has a voltage drop of Q/C. When switch S is closed at time t = 0, the sum of the
voltage across the resistor and the capacitor must equal the applied voltage. This gives us
the equation
1
RI + Q = E(t) .
C
The current in the circuit is the rate of change of the amount of charge on the capacitor. So
using the relationship dQ/dt = I, this becomes a first order differential equation for Q(t):
R

dQ
1
+ Q = E(t) .
dt
C
1

The initial condition for this equation is Q0 = Q(0), the initial amount of charge on the
capacitor. (Q0 could be set by imposing a voltage V0 = Q0 /C across the capacitor before
inserting it into the circuit.) In this lab, we will use Euler’s method to numerically solve this
differential equation for two different applied voltages: a constant voltage E0 , and then an
AC voltage E(t) = 117 sin(120πt) which corresponds to the voltage out of a standard wall
socket in the United States.

Prelab assignment
Before arriving in lab, answer the following questions. You will need your answers in lab to
work 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. Verify
that the function
Q1 (t) = E0 C 1 − e−t/(RC)
(1)
is a solution to the initial value problem with E(t) = E0 ,
1
1
dQ1
=−
Q1 + E0 ,
dt
RC
R

Q1 (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 that
Q1 (t) satisfies the differential equation. (In other words, do not try to find the
solution of the initial-value problem, but rather just check that the given function
solves 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” on
Table 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 function
Q2 (t) =

E0
120π e−γt − cos(120πt) + γ sin(120πt)
2 + (120π)2 )
R (γ

(3)

satisfies the initial-value problem with this E(t),
dQ2
1
1
=−
Q2 + E0 sin(120πt) ,
dt
RC
R

Q2 (0) = 0 ,

(4)

where the constant γ = 1/(RC) has units of Hz=sec−1 . What do the initial
1
condition and forcing term R E0 sin(120πt) represent for the system at the time
the switch is closed?

2

2. Suppose you implement Euler’s method using Matlab, using step size h, and create
a vector t of time steps from t = 0 to t = 1. Often we refer to the first entry as
t0 = 0, the next as t1 and the final entry will be tN = 1 where N h = 1. Matlab does
not enumerate these entries in the same way. The first element of a vector in Matlab
is always t(1). In this case, we will have t(1)=0, and t(N+1)=1. Find the Matlab
indices 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 = 2400
Record these values of n in Table 1 below.
N
12
120
1200
2400

n so that t(n)=0

n so that t(n)=.5

n so that t(n)=.86

n so that t(n)=1

Table 1: Matlab Indices for Time Vector

In the lab
The primary manner in which we will work with Matlab is by writing simple programs that
we can then run in Matlab. Each program is a small file with a .m extension. For this lab we
will write programs that implement Euler’s method, and will compare the results obtained
for several step sizes h with exact solutions (when these are available).
Creating a (new) file
Launch Matlab. Use the buttons in the “Current Folder” box to the left of the Matlab
window to make sure that you are working on the Desktop. Note that we will work with files
on 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 use
after 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 window
called untitled will open; give it a name by selecting File→Save As in this window; give your
file the name EULER.m and save it. You have now created an (empty) file called EULER.m.
Note that capitalization matters here.

3

The contents of the EULER.m file
This program will implement Euler’s method to solve the initial-value problem
dy
= f (t, y) ,
dt

y(a) = y0 .

(5)

In this application, the function y(t) stands for Q(t), the charge on the capacitor. In the
program below, everything following a % is a comment. Comments give you information
about the program, but are not evaluated by the computer. You may choose whether or not
to type these comments into your program, but if you include the comments in your file you
must include the %, or the computer will try to read them. Type the following program into
the file EULER.m you have just created:
clear t
% Clears old time steps and
clear y
% y values from previous runs
a=0;
% Initial time
b=1;
% Final time
N=12;
% Number of time steps
y0=0;
% Initial value y(a)
h=(b-a)/N;
% Time step
t(1)=a;
y(1)=y0;
for n=1:N
% For loop, sets next t,y values
t(n+1)=t(n)+h;
y(n+1)=y(n)+h*f(t(n),y(n)); % Calls the function f(t,y)=dy/dt
end
plot(t,y)
title([’Euler Method using N=’,num2str(N),’ steps, by MYNAME’])
% Include your own name
Substitute your own name in place of MYNAME so it will appear on the plots you create. When
you are finished typing in the program statements, save your work by choosing Save from
the 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.m
Since the program EULER.m refers to the function f(t,y), we must create a second .m file
called f.m to define the right-hand side of the differential equation for Matlab. Go to the File
menu 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. Replace
that with f.m and click Save as you did before. You should then type into the file f.m the
following commands:

4

function f=f(t,y)
R=18000;
% Resistance
C=.0000125;
% Capacitance
E0=117;
% Constant Voltage
f=-y/(R*C)+E0/R;
% Defines the function f
Be 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 differential
equation of the initial-value problem (5) above; the particular form of f (t, y) appearing on the
last line corresponds to the specific initial-value problem (2). To solve a different differential
equation with EULER.m or another solver, you need only change this file. The other lines in
f.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 the
case of a constant forcing E(t): equation (2). We will change this later to solve in the case
of a non-constant forcing.
The exact solution of the initial-value problem and the file yE.m
Since for this problem we also know a formula for the exact solution to the initial-value
problem, we can write a short Matlab program to evaluate this formula so we can compare
with the results of using Euler’s method for various step sizes h. To tell Matlab about the
exact solution formula, you should create another new file called yE.m. The file yE.m will
contain 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 equation
with EULER.m or one of the other numerical methods we will study in Lab 3, and you wish
to compare with an analytical expression for the exact solution, you should modify the file
yE.m as well as f.m. Type the following commands into a new file called yE.m, saving your
work when you have finished:
function yE=yE(t)
R=18000;
% Resistance
C=.0000125;
% Capacitance
E0=117;
% Constant Voltage
gam=1/(R*C);
yE=E0*C*(1-exp(-gam*t)); % Exact solution yE

Running your program
So far you have written a program and saved it. Now you want to use the program. In the
Command Window type EULER. The plot of the solution curve should then appear in a new
window.
5

Here is a summary of the programming we have just done. To solve an initial-value
problem 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 using
Euler’s method. Then you must tell your solver what it is supposed to solve; that is your
f.m file defining the right-hand side of the differential equation. Finally, the solver must be
told where to start, so you must specify initial conditions. In this case the initial conditions
are entered directly into the solver with the lines a=0 and y0=0. You ‘fine tune’ your solver
by changing N. The stopping time is also entered in the solver as b=1.

Lab problems
1. 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 on
the figure window and choose Print. . . ; a dialog box will appear. If you click OK
in 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 size
of your plot, add text, labels, legends, titles, etc.
(b) Using your indices from Prelab problem 2 and your approximate solution vector
y 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 type
y(n) in the Matlab command window and press return. Similarly, to view the
corresponding exact solution value, type yE(t(n)) and press return. Make sure
you have enough significant digits to compare the approximate and exact answers.
You can see more significant digits in Matlab by typing format long in the
Matlab command window and then hitting return.
2. Repeat Problem 1 again using N = 120 and N = 1200. Use these data to complete
the 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 of
a non-constant forcing function. To solve this problem using EULER.m, you will have
to change your file f.m so that the function value is different. Change the line
f=-y/(R*C)+E0/R;
to the line
f=-y/(R*C)+(E0/R)*sin(120*pi*t);
Note: probably the best way to do this is to comment the previous definition, but to
leave it in the file, so that your file f.m becomes

6

function 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);

%
%
%
%
%

Resistance
Capacitance
Constant Voltage
Old definition of f
New definition of f

To compute the exact solution Q2 (t) for the initial-value problem (4), you will similarly
need to modify the file yE.m. Change the line
yE=E0*C*(1-exp(-gam*t);
to
% yE=E0*C*(1-exp(-gam*t);
And then add after that the lines
omg=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 at
the bottom of the window.
To plot the approximate solution and the exact one on the same set of axes, replace
the line plot(t,y) near the end of the file EULER.m with the lines
plot(t,y,’:’, t,yE(t),’-’)
legend(’Approximate’,’Exact’)
This defines a vector corresponding to the exact solution. Sample values of the exact
solution are plotted connected by a solid line, and the Euler’s method approximation
is plotted on the same axes as a dotted line. The legend command puts a helpful
box in the corner of the plot to help you identify which graph corresponds to which
function. You can change the location of the legend on the graph by clicking on it or
by using the Edit menu on the figure.
5. Run EULER using N = 120, N = 1200, and N = 2400. Be sure that you understand
what is being graphed: for N = 120, we are plotting the Euler approximation to the
solution at the points t = 0, t = 0.008333, t = 0.016667, etc.—and are plotting the
value of the exact solution, yE , at the same points. Explain how this results in the
graph 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 the
accuracy of the numerical method, which we saw when we considered the constant
forcing function—as we increased the number of steps, the accuracy of the approximation improved. Another, undersampling, comes in to play when an insufficient number
of points are used to graph the solution. Undersampling occurs when the data points
are spaced too far apart to capture all the behavior of the equation. You may know
7

this phenomenon as “aliasing,” when an insufficiently sampled high frequency signal
appears to be a lower frequency signal. Even if you have calculated the exact solution
correctly its graph may not be accurate. Imagine plotting a sine wave with frequency
120 cycles per second by sampling the function 120 times a second; the sample points
will suggest a constant function instead of an oscillatory function. Did you see this
problem 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 prelab
problems, your solutions to lab problems 1–6, and a brief, original, conclusions
paragraph summarizing what you have learned.

Time

Index

Exact y

Approximate y

Error

t

n

yE(t)

y(n)

|y(n)−yE(t)|

Percent Error
y(n) − yE(t)
× 100
yE(t)

N = 12
.5
1
N = 120
.5
1
N = 1200
.5
.86
1
Table 2: Approximate Solutions to Equation (2)

last modified 9 Sep 2015

8

Get a 30 % discount on an order above $ 100
Use the following coupon code:
RESEARCH
Order your essay today and save 30% with the discount code: RESEARCHOrder Now
Positive SSL