The first guess is using the forward Euler approach. The first jacobian is using finite differences. Then NR method is used to solve for the next iteration and Broyden's method is used to update the jacobian until tolerance is met. The figure shows the results that I got compared with correct results. What improvements can I make to my code so that I am moving in the right direction(at the very least). Here is the MATLAB code.
function [t,y] = bdf1_qn_broyden(fcn,tspan,Y0,h,tol)
t = (tspan(1):h:tspan(2))'; %timestep vector
n = length(t); %no of steps
y = zeros(n,length(Y0)); %creting output array
y(1,:) = Y0; %first row of o/p is initial condition
i = 2;
jac = zeros(length(Y0));
unity = eye(length(Y0));
while i <= n
H = h*unity;
for j = 1:1:length(Y0)
jac(:,j) = fcn(t(i),y(i-1,:)' + H(:,j)) - fcn(t(i),y(i-1,:)');
end
iterate = (unity - jac)\unity;
jac = (1/h)*jac;
ycurr = y(i-1,:)' + h*fcn(t(i-1),y(i-1,:)');
count = 0;
diff = 1;
while any(diff > tol) && (count < 25)
curreval = fcn(t(i),ycurr);
ynext = ycurr - iterate*(ycurr - y(i-1,:)' - h*curreval);
delta = ynext - ycurr;
deltaeval = fcn(t(i),ynext) - curreval;
jac = jac + (1/dot(delta,delta))*(deltaeval - jac*delta)*delta';
iterate = (unity - h*jac)\unity;
diff = abs(delta);
ycurr = ynext;
count = count +1;
end
if count >= 25
disp('Iterative method failed to converge within 25 steps');
end
y(i,:) = ynext';
i = i+1;
end
end
As you might've guessed, I am new to this and my code is quite inefficient. However, I would like to atleast be accurate and slow if I'm going to be slow. Is this not possible without use of adaptive step size? Is BDF1 inept at solving such problems? Should I choose a different method? Or should I solve my implicit function better? Is there a mistake in my implementation? I don't know which way to proceed. Please point me in the right direction. Huge Thanks!