#lab 6 solving ODEs IVPs


%we solve a simple scalar equations:
%y'=ay y(0)=1
global a
function y=f(x,t)
  global a
  y=a*x;
endfunction

a=1;
t=linspace(0,1,20);
X=lsode(@f,x0=1,t)
plot(t,X);pause(1);
norm(X-exp(a*t'),'inf')


a=-1;
t=linspace(0,1,20);
X=lsode(@f,x0=1,t)
plot(t,X);pause(1);
norm(X-exp(a*t'),'inf')

%solving backward
%try the simples:  just  take time backward
t=linspace(-1,0); tb=fliplr(t);
Xb=lsode(@f,1,tb);
X=flipud(Xb);
plot(t,X)
%works! in older octave's version it did not... 
%one had to change time variable to -t in ODE and 
%solve the new ODE and get back the solution 

%blow up eq on [0,2]
t=linspace(0,2);
X=lsode(@(x,t) x*x,1,t)
%fails... try on [0,1-eps]... eps smaller and smaller
%you have to run the code furher manually 

%linear penduulum y''=-y y(0)=1 y'(0)=0
function y=lpen(x,t)
  y=zeros(2,1);
  x
  y(1)=x(2);
  y(2)=-x(1);
end
t=linspace(0,2*pi);
X=lsode(@lpen,[1;0],t);
plot(t,X);pause(1);  
#we can plot trajectory: x(1) vs x(2)
plot(X(:,1),X(:,2));pause(1);


%nonlinear penduulum y''=-sin(y) y(0)=1 y'(0)=0
function y=nlpen(x,t)
  y=zeros(2,1);
  x
  y(1)=x(2);
  y(2)=-sin(x(1));
end
t=linspace(0,2*pi);
X=lsode(@nlpen,[1;0],t);
plot(t,X);pause(1);  
#we can plot trajectory: x(1) vs x(2)
plot(X(:,1),X(:,2));pause(1);

#stiff equations
A=[998 1998;-999 -1999]
eig(A)
lsode_options("integration method","stiff");
tic;lsode(@(x,t) A*x,[1;1],[0,10]);t(1)=toc
tic;lsode(@(x,t) A*x,[1;1],[0,100]);t(2)=toc
lsode_options("integration method","non-stiff");
tic;lsode(@(x,t) A*x,[1;1],[0,10]);tn(1)=toc
tic;lsode(@(x,t) A*x,[1;1],[0,100]);tn(2)=toc
%non-stiff solver (in lsode Adams-Bashford failed...)
tn(1)/t(1) % stiff solver was 41 times faster on [0,10]
clear all

%problem two artifical numerical blow-up
X=lsode(f=@(x,t) 200*x,x0=1,[0,10])
%lsode() failed... due to numerical overflow in double precision...

%we can change variables e.g. y=c*x i.e. scale the solution
%y'=c*x'=c*200*x=200*y but y(0)=c*x(0)=c
X=lsode(f,1e-12*x0,[0,10])
%it stopped at differet later time wht suggests an artifical blowup...

%back to a real blow-up x'=x*x x(0)=1 x(t)=-1/(t-1)
X=lsode(@(x,t) x*x,1 ,[0,2])
%y=cx y'=cx'=cx*x=y*y/c y(0)=cx(0)=c y(t)=-c/(t-1)
c=1e-9;X=lsode(@(y,t) y*y/c,c,[0,2])
%lsode fails but at wrong time moment it should before one...
