diff options
| author | MK2020 <mh4115@ic.ac.uk> | 2017-03-14 20:37:55 +0000 |
|---|---|---|
| committer | GitHub <noreply@github.com> | 2017-03-14 20:37:55 +0000 |
| commit | dbdb6d78af349048aba9a6c715fa74f4435f1c9b (patch) | |
| tree | 00129df57777f0048a0f23a7fd0c125d68c449e7 | |
| parent | f50531293d34479f3e7985f8ba0207ec8a0fc388 (diff) | |
| download | NumericalAnalysis-dbdb6d78af349048aba9a6c715fa74f4435f1c9b.tar.gz NumericalAnalysis-dbdb6d78af349048aba9a6c715fa74f4435f1c9b.zip | |
Add files via upload
27 files changed, 1031 insertions, 1 deletions
diff --git a/Part 1/Code that sucks!/Huen_script.m b/Part 1/Code that sucks!/Huen_script.m new file mode 100644 index 0000000..f9f4561 --- /dev/null +++ b/Part 1/Code that sucks!/Huen_script.m @@ -0,0 +1,41 @@ +function Huen_script (tf) %tf is the end time
+
+%initailise the circuits
+R = 0.5;
+L = 1.5*10^(-3);
+h = 0.00001; %step size
+
+
+%initailise the container
+
+N = round(tf/h); %number of iterations
+t = zeros(1, N);
+Vout = zeros(1, N);
+current = zeros(1,N);
+
+%input voltage
+% step function of 5 volt
+%Vin = @(t)5*heaviside(t);
+Vin = @(t)4*sin(2*pi*6000*t);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+plot(Vout);
+xlabel('T/s');
+ylabel('Vout/V');
+end
\ No newline at end of file diff --git a/Part 1/Code that sucks!/error_script.m b/Part 1/Code that sucks!/error_script.m new file mode 100644 index 0000000..8e24418 --- /dev/null +++ b/Part 1/Code that sucks!/error_script.m @@ -0,0 +1,101 @@ +%-----------------------------------------------------------------
+clear;
+ts = 0; % set initial value of x_0
+is = 0;
+h = 0.0001; % set step-size
+tf = 0.03; % stop here
+R = 0.5;
+L = 0.0015;
+
+A = 6;
+T = 0.00015;
+
+vin = @(t) A * cos(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+%numerical solution
+vout = vin(t) - iout * R;
+
+%obtaining the exact solution with favorite method
+i_Exact = (6/L)*((2*pi)/T * sin((2*pi)/T*t) + (R/L)*cos((2*pi)/T*t) - (R/L)*exp(-(R/L)*t))/((2*pi)/T^2 + (R/L)^2);
+exact = vin(t) - R*i_Exact;
+
+%error as a function of t
+error = exact - vout;
+
+figure(1);
+plot(t,error); %for midpoint
+
+figure(2);
+
+for ind=3:16 % choose these carefully
+ h1=2^(-ind); % set stepsize
+
+
+ [t, iout] = midpoint(func, ts, tf, is, h1);
+ exact = vin(t) - R*i_Exact(t) ;
+ vout = vin(t) - R*iout;
+ error = max(abs(exact - vout));
+
+ plot(log(h1),log(error),'b*'); % log/log plot stepsize vs error
+ xlabel({'Time', '(seconds)'});
+ ylabel({'Error'});
+ title('MIDPOINT');
+ hold on; % for next value to be plotted
+end
+
+
+
+%
+% %--------------------------------------------------------------------
+%
+% clear;
+% ts = 0; % set initial value of x_0
+% is = 0;
+% h = 0.0000005; % set step-size
+% tf = 0.03; % stop here
+% R = 0.5;
+% L = 0.0015;
+%
+% A = 6;
+% T = 0.00015;
+%
+% vin = @(t) A * cos(2*pi*t/T);
+% func = @(t, iout) (vin(t) - iout*R) / L; % define func
+% [t, iout ] = ralston(func, ts, tf, is, h);
+%
+% %numerical solution
+% vout = vin(t) - iout * R;
+%
+% %obtaining the exact solution with favorite method
+% i_Exact = @(t) (6/L)*((2*pi)/T * sin((2*pi)/T*t) + (R/L)*cos((2*pi)/T*t) - (R/L)*exp(-(R/L)*t))/((2*pi)/T^2 + (R/L)^2);
+% exact = vin(t) - R*i_Exact(t);
+%
+% %error as a function of t
+% error = abs(exact - vout);
+%
+% figure(3);
+% plot(t,error); %for ralston
+%
+%
+%
+% figure(4);
+%
+% for ind=3:16 % choose these carefully
+% h1=2^(-ind); % set stepsize
+%
+%
+% [t, iout] = ralston(func, ts, tf, is, h1);
+% exact = vin(t) - R*i_Exact(t) ;
+% vout = vin(t) - R*iout;
+% error = max(abs(exact - vout));
+%
+% plot(log(h1),log(error)); % log/log plot stepsize vs error
+% xlabel({'Time', '(seconds)'});
+% ylabel({'Error'});
+% title('V_{out} versus time (original function)');
+% hold on; % for next value to be plotted
+% end
+
+
diff --git a/Part 1/Code that sucks!/error_script_heun.m b/Part 1/Code that sucks!/error_script_heun.m new file mode 100644 index 0000000..8e2dff3 --- /dev/null +++ b/Part 1/Code that sucks!/error_script_heun.m @@ -0,0 +1,45 @@ +function error_script_heun (tf)
+
+ts = 0; % set initial value of x_0
+is = 0;
+h = 0.0001; % set step-size
+R = 0.5;
+L = 0.0015;
+
+
+A = 6;
+T = 0.00015;
+
+vin = @(t) A * cos(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun(func, ts, is, h);
+
+%numerical solution
+vout = vin(t) - iout * R;
+
+%obtaining the exact solution with favorite method
+i_Exact = @(t) (6/L)*((2*pi)/T * sin((2*pi)/T*t) + (R/L)*cos((2*pi)/T*t) - (R/L)*exp(-(R/L)*t))/((2*pi)/T^2 + (R/L)^2);
+exact = vin(t) - R*i_Exact(t);
+%error as a function of t
+error = abs(exact - vout);
+
+figure(1);
+plot(t,error); %for heun
+%loglog(error);
+
+figure(5);
+for ind=3:16 % choose these carefully
+ h1=2^(-ind); % set stepsize
+
+
+ [t, iout] = heun(func, ts, is, h1);
+ exact = vin(t) - R*i_Exact(t) ;
+ vout = vin(t) - R*iout;
+ error = max(abs(exact - vout));
+
+ plot(log(h1),log(error),'b*'); % log/log plot stepsize vs error
+ hold on; % for next value to be plotted
+end
+
+end
+
diff --git a/Part 1/Code that sucks!/general_script.m b/Part 1/Code that sucks!/general_script.m new file mode 100644 index 0000000..4b0ca78 --- /dev/null +++ b/Part 1/Code that sucks!/general_script.m @@ -0,0 +1,104 @@ +clear;
+ts = 0; % set initial value of x_0
+is = 0;
+h = 0.0001; % set step-size
+tf = 0.03; % stop here
+R = 0.5;
+L = 0.0015;
+
+vin = @(t) 3.5;
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun_temp(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (original function)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 3.5;
+tau = 0.000150;
+
+vin = @(t) A * exp(-t.^2/tau);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun_temp(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (exponential square funtion)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 3.5;
+tau = 0.000150;
+
+vin = @(t) A * exp(-t/tau);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun_temp(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (exponential function)');
+
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 4;
+T = 0.0015;
+
+vin = @(t) A * sin(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun_temp(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (sine function)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 4;
+T = 0.0015;
+
+vin = @(t) A * square(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun_temp(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (square function)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 4;
+T = 0.0015;
+
+vin = @(t) A * sawtooth(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = heun_temp(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (sawtooth function)');
+
+
+
+
+
+
diff --git a/Part 1/Code that sucks!/heun.m b/Part 1/Code that sucks!/heun.m new file mode 100644 index 0000000..aa1b2e7 --- /dev/null +++ b/Part 1/Code that sucks!/heun.m @@ -0,0 +1,11 @@ +function [x,y] = heun(func, xa, ya, h)
+
+%caluclaye next x buy incrementing by the stepsize
+x = xa + h;
+
+gradient1 = feval(func, xa, ya); %calculate the gradient at t
+ypredictor=ya+h*gradient1; %calculate predictor for the next value of y
+gradient2=feval(func, x, ypredictor); %calculate the gradient at t + h
+y = ya + h/2*(gradient1 + gradient2);
+
+end
\ No newline at end of file diff --git a/Part 1/Code that sucks!/zihan_error.m b/Part 1/Code that sucks!/zihan_error.m new file mode 100644 index 0000000..cbabc50 --- /dev/null +++ b/Part 1/Code that sucks!/zihan_error.m @@ -0,0 +1,37 @@ +clear;
+%inital
+ts = 0; % set initial value of x_0
+is = 0;
+
+%parameter
+tf = 0.03; % stop here
+R = 0.5;
+L = 0.0015;
+
+
+
+%vin and the function
+A = 6;
+T = 0.00015;
+vin = @(t) A * cos(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+
+
+for ind=3:16 % choose these carefully
+ h1=2^(-ind); % set stepsize
+
+ [t, iout, exact] = midpoint(func, ts, tf, is, h1);
+
+ figure;
+ plot(exact)
+
+% vout = vin - R*iout;
+%
+% error = max(abs(exact - vout));
+%
+% plot(log(h1),log(error),'b*'); % log/log plot stepsize vs error
+% xlabel({'Time', '(seconds)'});
+% ylabel({'Error'});
+% title('MIDPOINT');
+% hold on; % for next value to be plotted
+end
diff --git a/Part 1/Code that sucks!/zihan_error2.m b/Part 1/Code that sucks!/zihan_error2.m new file mode 100644 index 0000000..1e7956d --- /dev/null +++ b/Part 1/Code that sucks!/zihan_error2.m @@ -0,0 +1,51 @@ +clear;
+ts = 0; % set initial value of x_0
+is = 0;
+tf = 0.03; % stop here
+R = 0.5;
+L = 0.00015;
+
+A = 6;
+T = 0.00015;
+
+vin = @(t) A * cos(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+
+
+
+figure(1);
+
+for ind=16:20 % choose these carefully
+
+ h1=2^(-ind); % set stepsize
+
+ %obtaining the exact solution with favorite method
+ [t, i_Exact] = exact_solution(R,L,T,tf,h1);
+ exact = vin(t) - R*i_Exact;
+
+ [t, iout] = ralston(func, ts, tf, is, h1);
+ vout = vin(t) - R*iout;
+
+ %error as a function of t
+ error1 = exact - vout;
+ error_final = max(abs(error1));
+
+ plot(log(h1),log(error_final), '*g'); % log/log plot stepsize vs error
+ xlabel({'Step size', '(h)'});
+ ylabel({'Error', '(volts)' });
+ title('RALSTON ERROR (LOGLOG PLOT) ');
+ hold on;
+
+
+end
+
+%this gives us the error function for the above method
+
+figure(2);
+ plot(error1);
+ xlabel({'Time', '(sec)'});
+ ylabel({'Error', '(volts)'});
+ title('RALSTON ERROR');
+
+
+
\ No newline at end of file diff --git a/Part 1/Good Code/Figures_Ex2/Midpoint_error.png b/Part 1/Good Code/Figures_Ex2/Midpoint_error.png Binary files differnew file mode 100644 index 0000000..f71fe37 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/Midpoint_error.png diff --git a/Part 1/Good Code/Figures_Ex2/Ralston_error.png b/Part 1/Good Code/Figures_Ex2/Ralston_error.png Binary files differnew file mode 100644 index 0000000..99e40c3 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/Ralston_error.png diff --git a/Part 1/Good Code/Figures_Ex2/Thumbs.db b/Part 1/Good Code/Figures_Ex2/Thumbs.db Binary files differnew file mode 100644 index 0000000..7bf8fe6 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/Thumbs.db diff --git a/Part 1/Good Code/Figures_Ex2/heun_error.png b/Part 1/Good Code/Figures_Ex2/heun_error.png Binary files differnew file mode 100644 index 0000000..ceb5b61 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/heun_error.png diff --git a/Part 1/Good Code/Figures_Ex2/heun_error_loglog.png b/Part 1/Good Code/Figures_Ex2/heun_error_loglog.png Binary files differnew file mode 100644 index 0000000..076c8c5 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/heun_error_loglog.png diff --git a/Part 1/Good Code/Figures_Ex2/heun_midpoint_ralston_comparsion_error.png b/Part 1/Good Code/Figures_Ex2/heun_midpoint_ralston_comparsion_error.png Binary files differnew file mode 100644 index 0000000..01b7a53 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/heun_midpoint_ralston_comparsion_error.png diff --git a/Part 1/Good Code/Figures_Ex2/midpoint_error_loglog.png b/Part 1/Good Code/Figures_Ex2/midpoint_error_loglog.png Binary files differnew file mode 100644 index 0000000..8174556 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/midpoint_error_loglog.png diff --git a/Part 1/Good Code/Figures_Ex2/ralston_error_loglog.png b/Part 1/Good Code/Figures_Ex2/ralston_error_loglog.png Binary files differnew file mode 100644 index 0000000..5ad9f35 --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/ralston_error_loglog.png diff --git a/Part 1/Good Code/Figures_Ex2/ralston_midpoint_comparison_error.png b/Part 1/Good Code/Figures_Ex2/ralston_midpoint_comparison_error.png Binary files differnew file mode 100644 index 0000000..ac2c7aa --- /dev/null +++ b/Part 1/Good Code/Figures_Ex2/ralston_midpoint_comparison_error.png diff --git a/Part 1/Good Code/comparison_error.m b/Part 1/Good Code/comparison_error.m new file mode 100644 index 0000000..2fe5823 --- /dev/null +++ b/Part 1/Good Code/comparison_error.m @@ -0,0 +1,61 @@ +clear;
+%intial values
+ts = 0;
+is = 0;
+tf = 0.03;
+
+%params
+R = 0.5;
+L = 0.0015;
+
+
+%input vin
+A = 6;
+T = 0.00015;
+
+vin = @(t) A * cos(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+
+
+figure;
+for ind=16:28 % choose these carefully
+
+ h1=2^(-ind); % set stepsize
+
+ %obtaining the exact solution with favorite method
+ [t, i_Exact] = exact_solution(R,L,T,tf,h1);
+ exact = vin(t) - R*i_Exact;
+
+ [t, iout_ralston] = ralston(func, ts, tf, is, h1);
+ vout_ralston = vin(t) - R*iout_ralston;
+
+ error_ralston = exact - vout_ralston;
+ error_ralston = max(abs(error_ralston));
+
+ [t, iout_heun] = heun_temp(func, ts, tf, is, h1);
+ vout_heun = vin(t) - R*iout_heun;
+
+ error_heun = exact - vout_heun;
+ error_heun = max(abs(error_heun));
+
+ [t, iout_midpoint] = midpoint(func, ts, tf, is, h1);
+ vout_midpoint = vin(t) - R*iout_midpoint;
+
+ error_midpoint = exact - vout_midpoint;
+ error_midpoint = max(abs(error_midpoint));
+
+
+ % comparision graphs (Heun, Midpoint and Ralston)
+ hold on;
+ plot(log(h1), log(error_ralston), 'r*','DisplayName', 'Ralston(loglog)error');
+ plot(log(h1), log(error_heun), 'g*', 'DisplayName', 'Heun(loglog)error');
+ plot(log(h1), log(error_midpoint), 'b*','DisplayName', ' Midpoint(loglog)eror');
+
+ legend('show')
+
+ xlabel({'Step size', '(h)'});
+ ylabel({'Error', '(volts)'});
+ title('ERROR COMPARISON');
+ hold on;
+
+end
diff --git a/Part 1/Good Code/error_script.m b/Part 1/Good Code/error_script.m new file mode 100644 index 0000000..eac74dd --- /dev/null +++ b/Part 1/Good Code/error_script.m @@ -0,0 +1,57 @@ +clear;
+ts = 0; % set initial value of x_0
+is = 0;
+tf = 0.03; % stop here
+R = 0.5;
+L = 0.0015;
+
+A = 6;
+T = 0.00015;
+
+vin = @(t) A * cos(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+
+
+
+figure(1);
+
+for ind=16:25 % This increments the time step by 2^-n every time where n = 16, 17, .... 25
+ %we started from 16 as we observed the results with larger step sizes
+ %gave irregualar data
+
+ h1=2^(-ind); % set stepsize
+
+ %obtaining the exact solution with favorite method
+ [t, i_Exact] = exact_solution(R,L,T,tf,h1);
+ exact = vin(t) - R*i_Exact;
+
+ %call appropriate method
+ %[t, iout] = heun_methodforerror(func, ts, tf, is, h1);
+ %[t, iout] = ralston(func, ts, tf, is, h1);
+ [t, iout] = midpoint(func, ts, tf, is, h1);
+
+ vout = vin(t) - R*iout;
+
+ %error as a function of t
+ error1 = exact - vout;
+ error_final = max(abs(error1));
+
+ %this plots a log lgo plot to show the order of the error
+ plot(log(h1),log(error_final), '*g'); % log/log plot stepsize vs error
+ xlabel({'Step size', '(h)'});
+ ylabel({'Error', '(volts)' });
+ title('RALSTON ERROR (LOGLOG PLOT) ');
+ hold on;
+
+
+end
+
+%this gives us the error function for the above method
+figure(2);
+ plot(error1);
+ xlabel({'Time', '(sec)'});
+ ylabel({'Error', '(volts)'});
+ title('RALSTON ERROR');
+
+
+
\ No newline at end of file diff --git a/Part 1/Good Code/exact_solution.m b/Part 1/Good Code/exact_solution.m new file mode 100644 index 0000000..b55e0b5 --- /dev/null +++ b/Part 1/Good Code/exact_solution.m @@ -0,0 +1,13 @@ +%this is our exact solution using our favorite method into a function
+function [ t, i_Exact ] = exact_solution(R,L,T,tfinal, h)
+
+%initial arrays
+N = round((tfinal) / h);
+t = zeros(1,N);
+i_Exact = zeros(1,N);
+
+%for loop calculating next exact current
+for i = 1 : N - 1
+ t(i+1) = t(i) + h;
+ i_Exact(i+1) = (6/L)*((2*pi)/T * sin((2*pi)/T*t(i)) + (R/L)*cos((2*pi)/T*t(i)) - (R/L)*exp(-(R/L)*t(i)))/((2*pi)/T^2 + (R/L)^2);
+end
diff --git a/Part 1/Good Code/heun.m b/Part 1/Good Code/heun.m new file mode 100644 index 0000000..aa1b2e7 --- /dev/null +++ b/Part 1/Good Code/heun.m @@ -0,0 +1,11 @@ +function [x,y] = heun(func, xa, ya, h)
+
+%caluclaye next x buy incrementing by the stepsize
+x = xa + h;
+
+gradient1 = feval(func, xa, ya); %calculate the gradient at t
+ypredictor=ya+h*gradient1; %calculate predictor for the next value of y
+gradient2=feval(func, x, ypredictor); %calculate the gradient at t + h
+y = ya + h/2*(gradient1 + gradient2);
+
+end
\ No newline at end of file diff --git a/Part 1/Good Code/heun_methodforerror.m b/Part 1/Good Code/heun_methodforerror.m new file mode 100644 index 0000000..ed38271 --- /dev/null +++ b/Part 1/Good Code/heun_methodforerror.m @@ -0,0 +1,13 @@ +function [y, t] = heun_methodforerror(f,a,b,ya,n)
+h = (b - a) / n;
+halfh = h / 2;
+y(1,:) = ya;
+t(1) = a;
+for i = 1 : n
+ t(i+1) = t(i) + h;
+ g = f(t(i),y(i,:));
+ z = y(i,:) + h * g;
+ y(i+1,:) = y(i,:) + halfh * ( g + f(t(i+1),z) );
+end;
+
+%func,t0,tf,i0,step_size
\ No newline at end of file diff --git a/Part 1/Good Code/heun_script.m b/Part 1/Good Code/heun_script.m new file mode 100644 index 0000000..34e8ac9 --- /dev/null +++ b/Part 1/Good Code/heun_script.m @@ -0,0 +1,223 @@ +function heun_script (tf) %tf is the end time
+
+%initailise the circuits
+R = 0.5;
+L = 0.0015;
+h = 0.0001; %step size
+
+
+%initailise the container
+
+N = round(tf/h); %number of iterations
+t = zeros(1, N);
+Vout = zeros(1, N);
+current = zeros(1,N);
+
+%input voltage
+% step function of 3.5 volt
+Vin = @(t)3.5*heaviside(t);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+figure(1);
+plot(Vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('(Heaviside) V_{out} versus time');
+
+figure(8);
+plot(t, current);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (Heaviside)');
+
+%---------------------------------------------------------------------------------------------
+
+%initailise the circuits information at the top
+
+%input voltage
+tau = 0.000150;
+A = 3.5;
+
+Vin = @(t) A * exp(-t.^2/tau);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+figure(2);
+plot(Vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (Exponential square function #1)');
+
+%-----------------------------------------------------------------------------------------------
+
+%initailise the container information at the top
+
+%input voltage
+tau = 0.000150;
+A = 3.5;
+
+Vin = @(t) A * exp(-t/tau);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+figure(3);
+plot(Vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (Exponential #2)');
+
+
+% %-------------------------------------------------------------------------------------------
+
+
+%initailise the circuits information at the top
+
+
+N = round(tf/h); %number of iterations
+t = zeros(1, N);
+Vout = zeros(1, N);
+current = zeros(1,N);
+
+%input voltage
+% step function of 5 volt
+% T= 0.00015, 0.000015, 0.0004, 0.0011
+T = 0.0011;
+Vin = @(t)4*sin(2*pi*t/T);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+figure(4);
+plot(Vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.0011)(Sine wave)');
+
+%-----------------------------------------------------------------------------------------
+
+%initailise the circuits information at the top
+
+%input voltage
+% step function of 5 volt
+A = 4;
+T = 0.0011;
+
+Vin = @(t) A * square(2*pi*t/T);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+figure(5);
+plot(Vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.0011) (Square wave)');
+
+%--------------------------------------------------------------------------------------
+%initailise the circuits information at the top
+
+%input voltage
+% step function of 5 volt
+A = 4;
+T = 0.0011;
+
+Vin = @(t) A * sawtooth(2*pi*t/T);
+
+%the initial condition
+t(1) = 0;
+current(1) = 0;
+
+
+%the equation
+func = @(t,current) (Vin(t)-R*(current))/L; %Function input for difference method
+
+
+%Huen
+for j = 1 : N-1
+ [t(j + 1),current(j + 1)] = heun(func, t(j), current(j), h);
+ Vout(j + 1) = Vin(t(j)) - R*current(j); %Create Vout array from Iout, R and Vin
+end
+
+%plot
+
+figure(6);
+plot(Vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.0011) (Sawtooth wave)');
+
+end
\ No newline at end of file diff --git a/Part 1/Good Code/midpoint.m b/Part 1/Good Code/midpoint.m new file mode 100644 index 0000000..000f103 --- /dev/null +++ b/Part 1/Good Code/midpoint.m @@ -0,0 +1,22 @@ +function [ x, y ] = midpoint( f, t0, tfinal, y0, h)
+
+%calculate number of steps
+N = round((tfinal - t0) / h);
+%initialise the array
+ya = zeros(1,N);
+ta = zeros(1,N);
+ya(1) = y0;
+ta(1) = t0;
+
+for i = 1 : N - 1
+ ta(i+1) = ta(i) + h;
+ halfstep = ta(i) + 1 * h / 2;
+ gradient1 = f(ta(i), ya(i));
+ ypredict = ya(i) + 0.5 * h * gradient1;
+ gradient2 = f(halfstep, ypredict);
+ ya(i+1) = ya(i) + h * gradient2;
+end
+x = ta;
+y = ya;
+
+
diff --git a/Part 1/Good Code/midpoint_script.m b/Part 1/Good Code/midpoint_script.m new file mode 100644 index 0000000..4543e48 --- /dev/null +++ b/Part 1/Good Code/midpoint_script.m @@ -0,0 +1,110 @@ +clear;
+ts = 0; % set initial value of x_0
+is = 0;
+h = 0.0001; % set step-size
+tf = 0.03; % stop here
+
+R = 0.5;
+L = 0.0015;
+
+vin = @(t) 3.5;
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+figure(1);
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (Heaviside)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 3.5;
+tau = 0.00015;
+
+vin = @(t) A * exp(-t.^2/tau);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+figure(2);
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (Exponential square function #1)');
+
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 3.5;
+tau = 0.00015;
+
+vin = @(t) A * exp(-t/tau);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+figure(3);
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (Exponential function #2)');
+
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 4;
+T = 0.0011;
+
+vin = @(t) A * sin(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+figure(4);
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.0011) (Sine wave)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+A = 4;
+T = 0.0011;
+
+vin = @(t) A * square(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+figure(5);
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.0011) (Square wave)');
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+A = 4;
+T = 0.0011;
+
+vin = @(t) A * sawtooth(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = midpoint(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+figure(6);
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.0011) (Sawtooth wave)');
+
+
+
+
+
+
diff --git a/Part 1/Good Code/ralston.m b/Part 1/Good Code/ralston.m new file mode 100644 index 0000000..4fc5a2e --- /dev/null +++ b/Part 1/Good Code/ralston.m @@ -0,0 +1,25 @@ +function [xa,ya] = ralston(func, t0, tf,i0 , h)
+%calculate number of steps
+N = round((tf) / h);
+%initialise the array
+xa = zeros(1,N);
+ya = zeros(1,N);
+xa(1)=t0;
+ya(1)=i0;
+
+%param
+a1=1/3;
+a2=2/3;
+p1=3/4;
+q11=3/4;
+
+for i=1: N - 1
+ xtemp=xa(i);
+ ytemp=ya(i);
+
+k1=func(xtemp,ytemp);
+k2=func(xtemp+p1*h,ytemp+q11*k1*h);
+
+ya(i+1)=ytemp+(a1*k1+a2*k2)*h;
+xa(i+1)=xtemp+h;
+end
diff --git a/Part 1/Good Code/ralston_script.m b/Part 1/Good Code/ralston_script.m new file mode 100644 index 0000000..4f4ab7a --- /dev/null +++ b/Part 1/Good Code/ralston_script.m @@ -0,0 +1,104 @@ +clear;
+ts = 0; % set initial value of x_0
+is = 0;
+h = 0.0001; % set step-size
+tf = 0.03; % stop here
+R = 0.5;
+L = 0.0015;
+
+% vin = @(t) 3.5;
+% func = @(t, iout) (vin(t) - iout*R) / L; % define func
+% [t, iout ] = ralston(func, ts, tf, is, h);
+%
+% vout = vin(t) - iout * R;
+% plot(t,vout);
+% xlabel({'Time', '(seconds)'});
+% ylabel({'V_{out}', '(volt)'});
+% title('V_{out} versus time (original function)');
+% %____________________________________________________________________
+% h = 0.0001;
+% tf = 0.03;
+% figure;
+% A = 3.5;
+% tau = 0.00015;
+%
+% vin = @(t) A * exp(-t.^2/tau);
+% func = @(t, iout) (vin(t) - iout*R) / L; % define func
+% [t, iout ] = ralston(func, ts, tf, is, h);
+%
+% vout = vin(t) - iout * R;
+% plot(t,vout);
+% xlabel({'Time', '(seconds)'});
+% ylabel({'V_{out}', '(volt)'});
+% title('V_{out} versus time (exponential square funtion)');
+% %____________________________________________________________________
+% h = 0.0001;
+% tf = 0.03;
+% figure;
+% A = 3.5;
+% tau = 0.00015;
+%
+% vin = @(t) A * exp(-t/tau);
+% func = @(t, iout) (vin(t) - iout*R) / L; % define func
+% [t, iout ] = ralston(func, ts, tf, is, h);
+%
+% vout = vin(t) - iout * R;
+% plot(t,vout);
+% xlabel({'Time', '(seconds)'});
+% ylabel({'V_{out}', '(volt)'});
+% title('V_{out} versus time (exponential function)');
+
+%____________________________________________________________________
+h = 0.0001;
+tf = 0.03;
+figure;
+A = 4;
+T = 0.000015;
+
+vin = @(t) A * sin(2*pi*t/T);
+func = @(t, iout) (vin(t) - iout*R) / L; % define func
+[t, iout ] = ralston(func, ts, tf, is, h);
+
+vout = vin(t) - iout * R;
+plot(t,vout);
+xlabel({'Time', '(seconds)'});
+ylabel({'V_{out}', '(volt)'});
+title('V_{out} versus time (T=0.000015) (sine function)');
+% % % ____________________________________________________________________
+% h = 0.0001;
+% tf = 0.03;
+% figure;
+% A = 4;
+% T = 0.0011;
+%
+% vin = @(t) A * square(2*pi*t/T);
+% func = @(t, iout) (vin(t) - iout*R) / L; % define func
+% [t, iout ] = ralston(func, ts, tf, is, h);
+%
+% vout = vin(t) - iout * R;
+% plot(t,vout);
+% xlabel({'Time', '(seconds)'});
+% ylabel({'V_{out}', '(volt)'});
+% title('V_{out} versus time (T=0.0011) (square function)');
+% %____________________________________________________________________
+% h = 0.0001;
+% tf = 0.03;
+% figure;
+% A = 4;
+% T = 0.000015;
+%
+% vin = @(t) A * sawtooth(2*pi*t/T);
+% func = @(t, iout) (vin(t) - iout*R) / L; % define func
+% [t, iout ] = ralston(func, ts, tf, is, h);
+%
+% vout = vin(t) - iout * R;
+% plot(t,vout);
+% xlabel({'Time', '(seconds)'});
+% ylabel({'V_{out}', '(volt)'});
+% title('V_{out} versus time (T=0.0011)(sawtooth function)');
+
+
+
+
+
+
@@ -1 +1,2 @@ -# Maths Numerical Analysis of ODEs using Matlab - Coursework 2017 +# Part 1 +Includes exercise 1 and 2 |
