Solved–Programming Assignment # 1– Solution

Description

Objective and Background

This assignment has two basic objectives:

1. To remind you that you know how to program (or that you should know!)

2. To give you a little practice solving di erential equations numerically

In this assignment you will write a program in C or C++ to solve rst-order di erential equations of the form

(D a)y₀(t) = 0; (1)

subject to some initial condition on y₀(t). One way (not the best, but probably the easiest) is to make the following approximation:

Dy0(t)	y0(t + t)		y0	(t)
					:	(2)
		t

Observe that in the limit as t ! 0, this approximation becomes exact. Substituting (2) into (1) we obtain

	y0(t + t) y0	(t)
			ay0	(t) = 0
	t
Solving for y0(t + t),
	y0(t + t) = (1 + a t)y0(t)				(3)

If we know y₀(t) at time t = 0, say y₀(0) = y₀ and we choose some t, then we can nd the (approximate)

value at later times by iterating (3):

With this de nition, it should be obvious that
x2(t) = x1(t)	(4)
x3(t) = x2(t)	(5)
and the entire di erential equation can be written as:
x3(t) + a2x3(t) + a1x2(t) + a0x1(t) = 0:	(6)

Let

• 3

x₁(t)

x(t) = 4 x₂(t) 5

x₃(t)

and let

• 3

^a11 ^a12 ^a13

A = 4 a₂₁ a₂₂ a₂₃ 5 :

^a31 ^a32 ^a33

Then a system of rst-order di erential equations

x(t) = Ax(t) (7)

may be written as

(D A)x(t) = 0

where the vector derivative is de ned as

• 3

x₁(t)

Dx(t) = 4 x₂(t) 5 :

x₃(t)

Then following the same steps as before, it is possible to approximate the next step of the solution:

x(t + t) = (I + A t)x(t): (8)

where I is the identity matrix,

• 3

1 0 0

I = 4	0	1	0	5 :
	0	0	1

How do we form the matrix A? Substituting (4), (5), and (6) into (7), we nd that A is given by

• 3

A = 4	0	1	0	5 :
	0	0	1
	a0	a1	a2

This structure can be extended to any order di erential equation.

Assignment

1. For the di erential equation

(D + 2:5)y₀(t) = 0

with initial condition y₀(0) = 3:

1. Find and plot the analytical (exact) solution to the di erential equation for 0 t 10.

2. Write a program in C(++) to plot a numerical solution using (3). You may have to try several values of t to get a good enough approximation.

2

1. 3. Compare the exact solution with the approximate solution.

2. For the third-order di erential equation

(D³ + 0:6D² + 25:1125D + 2:5063)y₀(t) = 0

with initial conditions y₀(0) = 1:5, y₀(0) = 2, y(0) = 1:

1. 1. Find and plot the analytical solution to the di erential equation for 0 t 10. Identify the roots of the characteristic equation and plot them in the complex plane.

1. 2. Put the third-order di erential equation into state-space form.

1. 3. Write a program in C(++) to plot an approximate solution using (8). You may have to try several values of t to get a reasonable approximation.

1. 4. Compare the exact solution with the approximate solution.

3. For the circuit shown here:

• ^R2

i1	(t)		+
		C	y(t)
f (t) +
		R1

where R₁ = 1 k , R₂ = 22 k , C = 10 F , and L = 5 H.

1. Determine a di erential equation relating the input f(t) to the output y(t).

2. Determine the initial conditions on y(t) if i₁(0) = 0:2 A and y(0) = 5 V.

3. Determine the analytical solution for the zero-input response of the system with these initial conditions.

4. Represent the di erential equation for the circuit in state variable form.

5. Using your program, determine a numerical solution to the di erential equation for the zero-input response.

6. Plot and compare the analytical and the numerical solution. Comment on your results.

7. Suppose that the circuit had nonlinear element in it, such as dependent sources. Describe how the analytical solution and numerical solution would change.

Turn in your programs, your plots, and your observations.

Hints and helps

Making plots To make plots, I recommend you write out an ascii le with the t and y values that you want. Then use Matlab to make the plots.

To write the les using C++, at the top of your program you need to open a le to write into:

#include <fstream> /* appears near the beginning of your program */

.

.

.

main()

{

.

.

.

ofstream outfile(“junk”); /* open a filed called ’junk’ for writing */

.

.

.

Then any time you want to write, you simply use the ofstream operators:

outfile << t << ” ” << y << endl;

At the end of the program you should close the le:

outfile.close();

The steps above will save the computed data into a le named \junk”. Then you can use the plot facilities of Matlab to make the plot. This is done by reading the data into Matlab, then plotting the data. The following commands are to be done inside Matlab.

load junk	% load the file	’junk’	into a variable called ’junk’
plot(junk(:,1),junk(:,2)) %		plot the first column vs. the second
	%	column

Matrix/vector operations The following are some simple, xed size matrix and vector routines that might be helpful.

void mattimes(double t, double a[][3],double m[][3])

{

int i,j;

for(i = 0; i < 3; i++)

for(j = 0; j < 3; j++)

m[i][j] = t*a[i][j];

}

void matsub(double a[][3],double b[][3],double c[][3])

{

int i,j;

for(i = 0; i < 3; i++)

for(j = 0; j < 3; j++)

c[i][j] = a[i][j]-b[i][j];

}

void matvecmult(double m[][3], double *v, double *prod)

{

double sum;

int i,j;

for(i = 0; i < 3; i++) {

sum = 0;

for(j = 0; j < 3; j++) {

sum += m[i][j]*v[j];

}

prod[i] = sum;

}

}

And here is a suggested way to use some of these (but check the signs carefully!):

double a[3][3] = {{0,-1,0},{0,0,-1},{20.002,100.05,0.4}}; double m[3][3];

double I[3][3] = {{1,0,0},{0,1,0},{0,0,1}};

…

mattimes(deltat,a,m);

matsub(I,a,m);

Miscellaneous hints Remember that the starting index for arrays in C/C++ is 0.

An example loop for the time variable:

for(t = 0; t <= 4; t += deltat) {

5