clear all;
close all;
Delta_t=0.1;
t=0:Delta_t:50;
%
%Preallocating Memory
%
x=zeros(length(t),3);
x(1,1)=0.0;x(1,2)=1.0;x(1,3)=1.0;
for n=1:length(t)-1
x(n+1,:)=Findxnplus1(Delta_t,x(n,:));
end
%Check with Matlab solution
[tmat,xmat] = ode23(@Example22p2,[0 50],[0 1 1]);
hold off
plot(t,x(:,1),’k-.’,tmat,xmat(:,1),’k-‘);
hold on
plot(t,x(:,2),’b-.’,tmat,xmat(:,2),’b-‘);
legend(‘x_0(t) Implicit Euler’,’x_0(t) MATLAB ode23()’,’x_1(t) Implicit Euler’,’x_1(t) MATLAB ode23()’);
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop at least once
p=xn; %Guess value for p
while eps>1.0e-13
[gfunc,jacobian]=Example22p2WithJacobian(p,Delta_t,xn);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient than inv(jacobian)*(-gfunc)
p=temp’+p;
end
end
function [gfunc,jac]=Example22p2WithJacobian(p, Delta_t,xn)
gfunc=[p(1)-xn(1)-p(2)*p(3)*Delta_t; p(2)-xn(2)+p(1)*p(3)*Delta_t; p(3)-xn(3)+0.5*p(1)*p(2)*Delta_t];
jac=[1 -p(3)*Delta_t -p(2)*Delta_t; p(3)*Delta_t 1 p(1)*Delta_t;0.5*p(2)*Delta_t 0.5*p(1)*Delta_t 1];
end
function dxdt = Example22p2(t,x)
dxdt = [x(2)*x(3); -x(1)*x(3);-0.5*x(1)*x(2)];
end