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

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

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 ﬁ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

to type these comments into your program, but if you include the comments in your ﬁle you

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’])

% Include your own name

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

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 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

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 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

**30 %**discount on an order above

**$ 100**

Use the following coupon code:

RESEARCH