Solving stiff ode system with vector parameters.

Hello ,

I'm trying to solve this stiff ode system , but the issue is that as long as the 'Ta' and 'qw' parameters are set as constants the code runs perfectly and returns good results, but when I change its values to be set as vectors , the cod starts to diverge and returns the following erreur, (Warning: Failure at t=7.190000e+02.Unable to meet integration tolerances without reducing the stepsize below the smallest value allowed) , Note that , the failure is at t=719 and its the same value of N-1 where N represents the length of the data vectors Ta0 and qw0), i have allready tried ode23s and ode15s, and I don't know any more tricks to overcom that so and help you can give is greatly appreciated.

%Solve With Ta qw set as constants. cod is running perfectly
clc;clear all
load life.dat
Ta0=life(:,7);  
qw0=life(:,4);

p=numel(Ta0);
x = 0:p-1;   
% Ta=@(t) interp1(x(:), Ta0(:), t); 
% qw=@(t) interp1(x(:), qw0(:), t); 
Ta=20;          %here I changed the values of Ta and qw with constants       
qw=350;

%step size
h=0.001;
N=ceil(p/h);
tspan=[0 N-1];

%parameters
q=250;
hi=4;
s=2.8*3;
l=0.1;
k=1.5;
fc=1;
c=920*850*s*l*fc;
Gc=hi*s;
Gd=k*s/l;
Gr=20;
Ti=25;
%the system of differential equations

   Eqns=@(t,T) ...
            [(Gc*(Ta-T(1))+Gr*(Ta-T(1))+qw*s)/c...
          ;(Gd*(T(1)-T(2))+Gc*(Ti-T(2)))/c];
%initial conditions           
T(:,1)=[20,20]; 
opts = odeset('RelTol',5.421011e-6,'AbsTol',5.421011e-6);
tic,[T0,X]=ode23s(@(t,T)Eqns(t,T), tspan, T,opts);toc
%plot the solution
figure
time1=linspace(datetime(2020,06,1,0,0,0),datetime(2020,07,01,0,0,0),length(T0));
plot(time1, X)
grid

the second cod:

%Solve With Ta qw set as vectors. cod is running but still diverges
clc;clear all
load life.dat
Ta0=life(:,7);  %varies from 10 to 30 
qw0=life(:,4);   %qw0 makes a sinusoidal fucntion and varies from 0 to 800 

%interpolation of Ta and qw 
p=numel(Ta0);
x = 0:p-1;   
Ta=@(t) interp1(x(:), Ta0(:), t); 
qw=@(t) interp1(x(:), qw0(:), t); 

%step size
h=0.001;
N=ceil(p/h);
tspan=[0 N-1];

%parameters
q=250;
hi=4;
s=2.8*3;
l=0.1;
k=1.5;
fc=1;
c=920*850*s*l*fc;
Gc=hi*s;
Gd=k*s/l;
Gr=20;
Ti=25;

%the system of differential equations
   Eqns=@(t,T) ...
            [(Gc*(Ta(t)-T(1))+Gr*(Ta(t)-T(1))+qw(t)*s)/c...
          ;(Gd*(T(1)-T(2))+Gc*(Ti-T(2)))/c]; 
 %initial conditions           
T(:,1)=[20,20]; 

opts = odeset('RelTol',5.421011e-6,'AbsTol',5.421011e-6);
tic,[T0,X]=ode15s(@(t,T)Eqns(t,T), tspan, T,opts);toc
%plot the solution
figure
time1=linspace(datetime(2020,06,1,0,0,0),datetime(2020,07,01,0,0,0),length(T0));

plot(time1, X)
grid

Answer 1

[darova]

I found the problem

[malek chek]

please explain further, if you are mentioning that i need to reduce the time step 'h' to 1, I'm afraid this doesn't work , because it also returns a divergence.

[darova]

tspan should be inbetween [0 .. p-1]

Can you attach the data so i can run the script and see what happens?

[malek chek]

Yes i will be grateful if you do that, i changed the file format to .txt so i can apload it , all i can help you with is the main characteristics that make it a very bad stiff equations is the data vectors, because they have sinusoidal format which varies strongly in compair to other parameters, finally the solution is supposed to have a simillar format to Ta0, and sorry for the dely .

[darova]

Here is the result i get

Is there something wrong with it?

[malek chek]

Yes, as you can see the solution that represents the first equation is diverging with time (the blue line) and as a result the other solution of the second equation (the grin line)is getting worse , but in reality the result should be something similar to the plot of Ta0, or you can also a hint on what the solution is supposed to be like if you reduce the parameter fc to 0.001 , but i'm i'm afraid this can be done just for illustration purposes ,

[darova]

What about constants?

c=920*850*s*l*fc /1000

[malek chek]

i'm afraid not , the value of c should be equal to 920*850*s*l*fc with fc set to 1,

[malek chek]

it's unusual type of stiffness differential equations, and needs a special solver and special care of some type to solve it. thank you again,