Solving stiff ode system with vector parameters.

Question

0 个投票

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