Attempting to use the trapezoidal rule to form a difference equation representing a circuit

I have a differential equation that has been proven to be correct.

The transfer function obtained by Laplace domain analysis and Matlab freqs match up and all is well.

The problem is somewhere in my attempt to digitize the differential equation. I am required to do this using the trapezium rule to approximate the analog signal.

The goal is to get a digitized equation in the form of..

Vo^n+1 = alpha1(Vi^n+1) + alpha2(Vi^n) + alpha3(Vo^n)

Where each alpha value is made up of a different combination of T, R1, R2 and C all numerical parameters.

I am following some guidelines set by my lecturer but if the approach seems wrong please let me know.

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Current incomplete solution..  Then I am trying to convert that equation to the Z domain and therefore use the freqz tool in Matlab. Please see the working out below. Here is the Matlab plot that shows what has gone wrong. The reason for the multiple lines are varying values of capacitance for C. The blue lines represent the validated analogue freqs(laplace domain transfer function) frequency response.

The green lines show my attempt at using the freqz(z domain discrete transfer function) frequency response.

If it is of any further interest, the code to generate those plots can be found below as well.

Any help would be sincerely appreciated, clutching at straws and the occasional hair at the moment!

Cheers

MATLAB code:

%% Comparing the impulse amlitude response using two different methods
%
%  Basically checking if freqs matches freqz

clc; clear all;

% Circuit Parameters
R1 = 2.2e3;
C = 0.01e-6;
R2 = 10e3;

fs = 10000;   % Sampling frequency (Hz)
T = 1/fs;     % Timestep

% Coefficients
den = (T*R1 + T*R2 + 2*R1*R2*C);

N = 1000;

% Digital parameters
y = zeros(N,1);
y(1) = 1;
result = zeros(N,1);

% Analogue parameters
fa = 0:fs/N:fs-fs/N;
w = 2*pi*fa;
j =sqrt(-1);

lastValue = 0;

clf;
figure(1);

hold on;

for i=1:5

% Analogue Recalculations
b = 1/(R1*C);
a = (R1 + R2)/(R1*R2*C);
Ha = freqs(b,[1 a], w);

% Digital Recalculations
alpha1 = (T*R2)/den;
alpha2 = alpha1;
alpha3 = (2*R1*R2*C - T*R1 - T*R2)/den;

[Hd,fd] = freqz([alpha1 alpha1],[1 alpha3],N,fs);

for n=1:N-1

currentValue = alpha1*y(n+1) + alpha2*y(n) - alpha3*lastValue;
result(n) = currentValue;
lastValue = currentValue;
end

C = 1.5*C;

plot(fa,abs(Ha),fd,abs(Hd)),xlabel('Frequency'),ylabel('H(s)');

end

hold off;
• Why bother plotting the analog response between 5k and 10k? Nothing in the digitized version will match that --- your sampling rate is only 10kHz, so only valid frequencies of the digital implementation are 0 to 5kHz. Will try to work through your equations... but not today. :-) – Peter K. Nov 29 '13 at 17:55