Numerical Solutions to Differential Equations

With Applications to Electric Circuits

Introduction to Numerical Differentiation

· The derivative of a time varying function can be approximated by the slope between two adjacent points (or the change in y value of the function, divided by the change in time).

· Thus, the derivative of a function y(t) can be approximated as follows:

· There is a function in MATLAB called diff( ) that allows you to obtain the differences above easily.

· Consider an array y with elements

y=[y(1), y(2),…, y(n)]

Calling diff(x) returns an array

[y(2)-y(1), y(3)-y(2),…, y(n)-y(n-1)]

· Try this

y=[3,2,6,3,4,5,1]

dy=diff(y)

· Note that dy is one element shorter than y and these differences are not divided by. To get the approximate derivative use diff(y)/deltat (deltat must be defined as the correct time spacing between samples in y).

· Similarly, a second derivative can be approximated as follows:

This can be achieved in MATLAB by using the diff( ) function twice diff(diff(y))/(deltat^2).

Computing Derivatives

· Let's dig in and compute some first and second derivative approximations for a function known as a sigmoid. Create a new m-file and type in the following code to generate a sigmoid function:

close all; clear; clc;

t=linspace(-10,10 ,1000);

y=1./(1+exp(-t));

plot(t,y)

xlabel('t ')

ylabel('y ')

title('Sigmoid Function ')

· Now let’s compute a numerical derivative

deltat=t(2)-t(1); % all evenly spaced

dir=diff(y)/deltat;

plot(t,y, t(1:length(dir)),dir)

legend('Sigmoid ', 'Derivative ')

title('Sigmoid Function and its Derivative ')

· Is this derivative what you expected? Note that if the original sigmoid represents the position of some moving object over time, the derivative function tells you the velocity of the object at each time. The second derivative would tell you the acceleration of the object.

· OK, so let's try a second derivative. Try to guess what the second derivative should look like before executing the code below. It will be the derivative of the first derivative.

deltat=t(2)-t(1); % all evenly spaced

dir1=diff(y)/deltat;

dir2=diff(diff(y))/deltat^2;

plot(t,y, t(1:length(dir1)),dir1, t(1:length(dir2)),dir2)

legend('Sigmoid', '1^st Derivative' , '2nd Derivative',2)

title('Sigmoid Function and its Derivatives')

Solving the RC Circuit

· Consider a flashbulb in a typical camera. In order to get the brief but powerful flash with just a small battery, we need to use a capacitor. A capacitor is simply two metal plates separated by an insulator. Charge builds up on this device slowly from the battery and this charge can be unleashed quickly into the flashbulb. Hooking the battery directly to the flash bulb would yield a very disappointing glow (at best).

· Let’s say we have a 9V battery in our camera and we charge up the capacitor so that it now has 9V across its terminals. We can now remove the battery (it has done its job) and we can connect the capacitor to the flashbulb when we are ready. The flashbulb can be modeled as a resistor. Thus, the circuit shown above models our situation nicely.

· We are interested to know what the voltage as a function of time will be after we trigger the flash, by connecting the charged up capacitor to the flash bulb. The voltage function will tell us how much “kick” the flash is getting and for how long. If we can develop a mathematical model for this process, we can easily test the effect of using different sized capacitors or different resistance bulbs. This will enable us to design a flash system to meet our needs.

· The capacitor and the resistor have the same voltage across them. Furthermore, they have the same magnitude of current flowing through them. What we need to do is figure out what voltage function v(t) can simultaneously exist across both components.

· The voltage-current relationship of a resistor is i(t)=v(t)/R. The voltage-current relationship for a capacitor in the configuration above is i(t)=-Cdv(t)/dt. Since these two components share i(t), the following differential equation holds true

Rearranging this, in a standard form for a differential equation, yields

· A function v(t) which satisfies this equation and starts at 9V (initial condition), will be the only v(t) that can exist in the circuit. We need to find a v(t) who’s derivative is –v(t)/(RC).

· Instead of trying to solve the difficult differential equation, let’s convert this into an equation with numerical approximations for the derivatives

· Solving for v(n) in terms of v(n-1) yields a recursive equation whereby we can solve for the next voltage point from the previous point. Knowing the process starts at 9 Volts is enough to solve for all future voltages

· Here is some MATLAB code to evaluate this recursive relationship

deltat=.0001; % in Seconds

C=.01; % Measured in Farads

R=.5; % Measured in Ohms

v(1)=9; % Volts

t(1)=0;

for n=2:1:1000

v(n)=( 1/(deltat/(R*C)+1) )*v(n-1);

t(n)=t(n-1)+deltat; % keep track of our time steps as well

end

plot(t,v)

xlabel('t (sec) ')

ylabel('v(t) (Volts) ');

title('RC Circuit Analysis ');

· Recall that power is voltage squared divided by resistance

figure

p= v.^2/R;

plot(t,p)

xlabel('t (sec) ')

ylabel('p(t) (Watts) ');

title('RC Circuit Analysis ');

· I hope that you are impressed that we can get over 150 Watts of power for our flash from a simple 9V battery! The catch is that we only get that power for a small fraction of a second.

· To see the numerical value of power at some particular time, you need to figure out which array element in p corresponds to the time you are interested in. Search the time array for the time you want. Let’s say you want the power at t=0.0020 sec. If you try looking at different elements in the t array, you will find that the 21’st element in t is 0.0020

t(21)

The power at time t=0.0020 sec. is then element 21 in p

p(21)

· To solve large scale differential equations with MATLAB see ode45 (i.e., type “help ode45” at the command window prompt).

· Here is an example: http://www.mit.edu/people/abbe/matlab/ode.html

· Here is more info on differential equations in MATLAB: http://www.mathworks.com/support/tech-notes/1500/1510.html