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/heun.m | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 coursework17/heun.m (limited to 'coursework17/heun.m') 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 + -- cgit v1.2.3-54-g00ecf