# Math 216: Diﬀerential Equations Lab 2: Euler’s Method and RC Circuits Goals

August 31, 2017

Question
Math 216: Diﬀerential 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 ﬁles 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 diﬀerential 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 ﬁrst order diﬀerential 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
diﬀerential equation for two diﬀerent 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
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 ﬁrst 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 satisﬁed, and then by diﬀerentiating the formula
(1) and comparing with the right-hand side of the diﬀerential equation show that
Q1 (t) satisﬁes the diﬀerential equation. (In other words, do not try to ﬁnd 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)

satisﬁes 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 ﬁrst entry as
t0 = 0, the next as t1 and the ﬁnal entry will be tN = 1 where N h = 1. Matlab does
not enumerate these entries in the same way. The ﬁrst 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 ﬁle 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) ﬁle
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 ﬁles
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 ﬁnish in the lab today. Anything you do not move to one of these will be lost.
Start editing a ﬁle 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
ﬁle the name EULER.m and save it. You have now created an (empty) ﬁle called EULER.m.
Note that capitalization matters here.

3

The contents of the EULER.m ﬁle
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
must include the %, or the computer will try to read them. Type the following program into
the ﬁle 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’])
Substitute your own name in place of MYNAME so it will appear on the plots you create. When
you are ﬁnished typing in the program statements, save your work by choosing Save from
the File menu, or by clicking on the ﬂoppy disk icon on the EULER.m editor window.
The right-hand side of the diﬀerential equation and the ﬁle f.m
Since the program EULER.m refers to the function f(t,y), we must create a second .m ﬁle
called f.m to deﬁne the right-hand side of the diﬀerential 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 ﬁle 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 ﬁnished typing.
The ﬁle f.m contains the function f (t, y) for the right-hand side of the general diﬀerential
equation of the initial-value problem (5) above; the particular form of f (t, y) appearing on the
last line corresponds to the speciﬁc initial-value problem (2). To solve a diﬀerent diﬀerential
equation with EULER.m or another solver, you need only change this ﬁle. The other lines in
f.m serve to deﬁne speciﬁc values for the constants; here we are using R = 18kΩ, C = 12.5µF,
and E0 = 117V.
Note that your newly created f.m ﬁle 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 ﬁle 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 ﬁle called yE.m. The ﬁle yE.m will
contain the exact solution Q1 (t) of the initial-value problem (2) corresponding to the righthand side function f (t, y) deﬁned in the ﬁle f.m. If you solve a diﬀerent diﬀerential 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 ﬁle
yE.m as well as f.m. Type the following commands into a new ﬁle called yE.m, saving your
work when you have ﬁnished:
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

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 diﬀerential equation you need a solver. In this case, EULER is your solver,
and generates a numerical approximation to the solution to the diﬀerential equation using
Euler’s method. Then you must tell your solver what it is supposed to solve; that is your
f.m ﬁle deﬁning the right-hand side of the diﬀerential 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 ‘ﬁne 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 ﬁgure 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 ﬁgure window to change the size
(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 signiﬁcant digits to compare the approximate and exact answers.
You can see more signiﬁcant 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 ﬁle f.m so that the function value is diﬀerent. 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 deﬁnition, but to
leave it in the ﬁle, so that your ﬁle 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 ﬁle 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 ﬁles 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 ﬁle EULER.m with the lines
plot(t,y,’:’, t,yE(t),’-’)
legend(’Approximate’,’Exact’)
This deﬁnes 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 ﬁgure.
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 diﬀerent when you use the diﬀerent values of N .
6. Many diﬀerent 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 insuﬃcient 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 insuﬃciently 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 modiﬁed 9 Sep 2015

8

Get a 30 % discount on an order above \$ 100
Use the following coupon code:
RESEARCH
Positive SSL