The backward Euler method for integrating the system x'=f(x) is
This is an implicit method: given x[i] you need to use a root-finding method, such as Newton's method, to find x[i+1]. Here is a matlab function bestep(x,h) that does exactly this. Put it in a file called bestep.m. The function f defining the system of differential equations should be in the file f.m. f and x should be taken to be column vectors.
function a=bestep(x,h) d=size(x,1); eps=0.00001; cha=eps*eye(d,d); M=zeros(d,d); fx=f(x); a=x+h*fx; con=1; while(con>0) fa=f(a); for i=1:d M(:,i)=(f(a+cha(:,i))-fa)/eps; end ch = (eye(d,d)-h*M)\(x+h*fa-a); con = norm(ch,1); a=a+ch; endNotes: The function uses Newton's method, with derivatives of f computed with the simplest finite difference formula with "small" parameter eps=0.00001. You might want to change this! It also continues to run the Newton iteration until it finds a y such that y=x+hf(y) to machine accuracy. You might want to change this too... just replace the statement "while(con>0)" with "while(con>1e-10)" or something similar.
Back to the course homepage
Back to my main teaching page
Back to my main page