blob: ab7edfd8acac2b4d5efc9737ee6039257c196e88 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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
|
