From abcd65d98b8a4830673d53eaa9665ad195b70b03 Mon Sep 17 00:00:00 2001 From: Vasil Zlatanov Date: Fri, 17 Feb 2017 14:43:50 +0000 Subject: work in progress solution for exercise 1 of coursework --- coursework17/error_script.m | 13 +++++++++++ coursework17/heun.m | 25 ++++++++++++++++++++++ coursework17/midpoint.m | 25 ++++++++++++++++++++++ coursework17/octave-workspace | 0 coursework17/ralston.m | 25 ++++++++++++++++++++++ coursework17/ralston_script.m | 50 +++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 138 insertions(+) create mode 100644 coursework17/error_script.m create mode 100644 coursework17/heun.m create mode 100644 coursework17/midpoint.m create mode 100644 coursework17/octave-workspace create mode 100644 coursework17/ralston.m create mode 100644 coursework17/ralston_script.m diff --git a/coursework17/error_script.m b/coursework17/error_script.m new file mode 100644 index 0000000..29988bf --- /dev/null +++ b/coursework17/error_script.m @@ -0,0 +1,13 @@ +% This script will carry out error analysis +T = 150e-6; +f = 1/T; + +A=6; + +data_points=10e3; + +% TODO find exact solution +exact = @(t) (); + +for j=1:data_points + my_error = exact(time_array(j)) - Vout_array(j); diff --git a/coursework17/heun.m b/coursework17/heun.m new file mode 100644 index 0000000..ab7edfd --- /dev/null +++ b/coursework17/heun.m @@ -0,0 +1,25 @@ +function [time_array, Vout_array] = heun(R, L, Vin, current_initial, step, t_final) + % Determine how many times to step forward + N = round(t_final/step); + + % Initialise some arrays for the output + current_array=zeros(1,N); + time_array=zeros(1,N); + current_array(1) = current_initial; + Vout_array(1) = Vin(0) - R*current_initial; + + % Given the input voltage as a function of time we can express + % dI/dt (dcurrent_dt) as a function of time as well (given R and L) + dcurrent_dt = @(t, I) (Vin(t) - I*R)/L; + + % Loop as many times as step allows and: + % calculate the gradient at point + % Use it to emtimate next two points + + for j=1:N-1 % loop N-1 times + grad = feval(dcurrent_dt, time_array(j), current_array(j)); + current_array(j+1)=current_array(j) + step/2*(grad + feval(dcurrent_dt, time_array(j)+step, current_array(j)+step*grad)) ; + time_array(j+1) = time_array(j) + step; + Vout_array(j+1) = Vin(time_array(j+1)) - R*current_array(j+1); + end + diff --git a/coursework17/midpoint.m b/coursework17/midpoint.m new file mode 100644 index 0000000..44ccdae --- /dev/null +++ b/coursework17/midpoint.m @@ -0,0 +1,25 @@ +function [time_array, Vout_array] = midpoint(R, L, Vin, current_initial, step, t_final) + % Determine how many times to step forward + N = round(t_final/step); + + % Initialise some arrays for the output + current_array=zeros(1,N); + time_array=zeros(1,N); + current_array(1) = current_initial; + Vout_array(1) = Vin(0) - R*current_initial; + + % Given the input voltage as a function of time we can express + % dI/dt (dcurrent_dt) as a function of time as well (given R and L) + dcurrent_dt = @(t, I) (Vin(t) - I*R)/L; + + % Loop as many times as step allows and: + % calculate the gradient at point + % Use it to emtimate next two points + + for j=1:N-1 % loop N-1 times + grad = feval(dcurrent_dt, time_array(j), current_array(j)); + current_array(j+1)=current_array(j) + step*feval(dcurrent_dt, time_array(j)+step/2, current_array(j)+step/2*grad); + time_array(j+1) = time_array(j) + step; + Vout_array(j+1) = Vin(time_array(j+1)) - R*current_array(j+1); + end + diff --git a/coursework17/octave-workspace b/coursework17/octave-workspace new file mode 100644 index 0000000..e69de29 diff --git a/coursework17/ralston.m b/coursework17/ralston.m new file mode 100644 index 0000000..d910442 --- /dev/null +++ b/coursework17/ralston.m @@ -0,0 +1,25 @@ +function [time_array, Vout_array] = ralston(R, L, Vin, current_initial, step, t_final) + % Determine how many times to step forward + N = round(t_final/step); + + % Initialise some arrays for the output + current_array=zeros(1,N); + time_array=zeros(1,N); + current_array(1) = current_initial; + Vout_array(1) = Vin(0) - R*current_initial; + + % Given the input voltage as a function of time we can express + % dI/dt (dcurrent_dt) as a function of time as well (given R and L) + dcurrent_dt = @(t, I) (Vin(t) - I*R)/L; + + % Loop as many times as step allows and: + % calculate the gradient at point + % Use it to emtimate next two points + + for j=1:N-1 % loop N-1 times + grad = feval(dcurrent_dt, time_array(j), current_array(j)); + current_array(j+1)=current_array(j) + step/4*(grad + 3*feval(dcurrent_dt, time_array(j)+2/3*step, current_array(j)+2/3*step*grad)) ; + time_array(j+1) = time_array(j) + step; + Vout_array(j+1) = Vin(time_array(j+1)) - R*current_array(j+1); + end + diff --git a/coursework17/ralston_script.m b/coursework17/ralston_script.m new file mode 100644 index 0000000..9db929f --- /dev/null +++ b/coursework17/ralston_script.m @@ -0,0 +1,50 @@ +% Vin = aI' + bI + +R = 0.5; % 0.5Ohm +L = 0.0015; % 1.5mH +data_points = 10000; + +% Go on for 8 time constants +time_constant = L/R; +step = time_constant*8/data_points; + +% Signal 1 Vin = 5.5V +Vin = @(t) 5.5 + 0*t; % 5.5V +current_initial=0; + +[t, Vout] = ralston(R, L, Vin, current_initial, step, data_points*step); +plot(t, Vout); + +% Signal 2a Vin = 5.5*exp(-t^2/r) V +Vin = @(t) 5.5*exp(-(t*160*10^(-6))^2) % 5.5V +current_initial=0; + +[t, Vout] = ralston(R, L, Vin, current_initial, step, data_points*step); +plot(t, Vout); + +% Signal 2b Vin = 5.5*exp(-t/r) +Vin = @(t) 5.5*exp(-t*160*10^(-6)) % 5.5V +current_initial=0; + +[t, Vout] = ralston(R, L, Vin, current_initial, step, data_points*step); +plot(t, Vout); + +% Signal 3 +T(1) = 20e-6; % 20us +T(2) = 450e-6; % 450us +T(3) = 1000e-6;% 1000us + +for j=1:3 + f = 1/T(j); + Vin = @(t) sin(2*pi*f*t); + [t, Vout] = ralston(R, L, Vin, current_initial, step, data_points*step); + plot(t, Vout); + + Vin = @(t) sawtooth(2*pi*f*t); + [t, Vout] = ralston(R, L, Vin, current_initial, step, data_points*step); + plot(t, Vout); + + Vin = @(t) square(2*pi*f*t); + [t, Vout] = ralston(R, L, Vin, current_initial, step, data_points*step); + plot(t, Vout); +end -- cgit v1.2.3