Why do I get different results from the numerical solution and the analytical solution with ODE15S? The numerical solution looks smoother than the analytical one

I want to simulate the dynamic behavior of pesticide in fruits, which is affected by daily temperature and relative humidity, by using an ode15s solver. And the result of the numerical model is then compared with the results of the analytical model.

Why do I get different results from the numerical solution and the analytical solution with ODE15S? The numerical solution looks smoother than the analytical one. And the results of ode15s solver seem to average the results to make it smoother, could somebody give me some ideas how to solve this problem? It gives me a headache. Thanks a lot!

[Alan Stevens]

Looks like the analytical values are calculated on a finer mesh than the numerical ones. However, it is difficult to tell without seeing both the analytical expression and the coding used to generate the numerical solution.

[Quanshun An]

Dear Alan Stevens,

The codes of the numerical model and the analytical model are attached, could you give me some idea how to solve this problem? Thanks!

Best regards,

Quanshun

[Alan Stevens]

Sorry! They are far too complicated for me to investigate, given that I don't know anything about the context. You also use a large number of global variables which never helps with debugging.

[Quanshun An]

Dear Alan,

The codes of the numierical model and the analytical model have been simplified and are attached; and I hope you can understand it.

Thanks a lot!

Best regard,

Quanshun

Answer 1

[Alan Stevens]

% replace
m_Soil = m_Soil_initial .* exp(- k1(i) .* t_measured(i+1));

% by
m_Soil = m_Soil_initial .* exp(- k1(i) .* (t_measured(i+1)-t_measured(i)));

[Quanshun An]

Dear Alan,

I really appreaicate your help, you are a genius!

Best regards,

Quanshun

Answer 2

[Alan Stevens]

Much easier to follow! Your "analytical" solution is incorrect for "constants" that change with time. See the following

%% Read data
[num, txt, raw] = xlsread('test','Sheet1');   % Read data from document
Time_Cell=raw(2:137,1);                                          % Extract time value
Time = cell2mat(Time_Cell);                                      % Convert cell to double 
Value_Cell=raw(2:137,2:3);                                       
Value = cell2mat(Value_Cell);                                    

t_measured = Time;                            
k1 = Value(:,1);                                        
b1 = Value(:,2);

m_Soil_initial = 10;
% Chemical mass
%m_Soil = b1 ./ k1 .* (1 - exp(- k1 .* t_measured)) + m_Soil_initial .* exp(- k1 .* t_measured);
m_Soil = m_Soil_initial .* exp(- k1 .* t_measured);
figure(1)
plot(t_measured, m_Soil,'k','LineStyle','-','LineWidth',1);
grid
axis([0 40 0 10])


% dm/dt = -k1*m    m = m0*exp(-k1*t)
% With k changing from one time to the next, you need to ingtegrate between
% the t_measured times starting from the previous time i.e. with a revised 
% "initial" value of mass.  See below 

% Integrate from 0 to 1
% m(t) = 10*exp(-k11*t)   m(1) = 10*exp(-k11)
% Integrate from 1 to 2
% m(t) = m(1)*exp(-k12*(t-1))   m(t) = 10*exp(-k11)*exp(-k12*(t-1)) 
%                               m(t) = 10*exp(-k11+k12)*exp(-k12*t)
%                               m(2) = 10*exp(-k11-k12)
% etc.
k = cumsum(k1);
m = [10; 10*exp(-k)]; 
m(end)=[];
hold on
plot(t_measured,m,'bo-')
legend('original','stepwise integration')

[Quanshun An]

Dear Alan,

Hope you are doing well!

Thanks for your great idea! I think you are right!

I modifed the code of analytical model according to your advice, but I got the different results, please see the code and figure.

Best regards,

Quanshun

clc; clear; close all; dbstop if error;

%% Read data
[num, txt, raw] = xlsread('test','Sheet1');   
Time_Cell=raw(2:137,1);                                          
Time = cell2mat(Time_Cell);                                      
Value_Cell=raw(2:137,2:3);                                       
Value = cell2mat(Value_Cell);                                    
t_measured = Time;                            
k1 = Value(:,1);                                        
b1 = Value(:,2);

% Calculation
m_Soil_initial = 10;
A = [];
for i = t_measured(2):t_measured(end)
    m_Soil = m_Soil_initial .* exp(- k1(i) .* t_measured(i+1));
    A = [A,m_Soil];
    m_Soil_initial = m_Soil;
end

A = [10,A]
figure(1)
plot(t_measured, A,'k','LineStyle','-','LineWidth',1);