fit 2nd order differential equation parameters to data

36 Ansichten (letzte 30 Tage)
I have a second order differential equation with some tunable parameters that I need to fit to some data. I have spent some time looking over some other postings with answers from (StarStrider especially, with some key assists from Torsten) and have finally gotten my script to run, but something isn't right. I suspect it's something I'm missing about passing in my y(t) data into the function that holds my system of differential equations, but I do not know for sure. I am hoping that someone can see what I'm missing. The data fitted function should follow the data much more closely than it does. The data is y vs t, and the diff. eq. is y''=a, where a is a function of y'. Here is my script:
clear all
close all
clc
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ]-.35 ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
a_m_guess = 200 ; %init guess for a_m, the parameter to be fitted
a_m = lsqcurvefit(@(am,t)FittingFunction(am,t),a_m_guess,Time_data,Y_data) ;
figure %plot the data as well as the fitted function
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(0,max(Time_data),100) ;
plot(plot_times, FittingFunction(a_m,plot_times),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest')
function F = FittingFunction(a_m,Time_data)
init_slope = 3.39 ; %This is a boundary condition
tspan = Time_data ;
InitConds = [0 init_slope] ; %y(t=0)=0, y'(t=0)=Us
[t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,a_m),tspan,InitConds) ;
F = y_atTdata(:,1) ; %return y only, not y and y'
end
function Y = SysOfDiffEqns(t,y,a_m)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs,
% substitute Y(1) = y
% Y(2) = y'
% makes a system of two first order ODEs:
% Y(1)' = Y(2)
% Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
D = 8.7049 ; % constant input parameter
U_m = 0.95*D ; % costant input parameter
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end
Thanks very much for any assistance. Again, I'm thinking I somehow need to pass in my y_data into @SysOfDiffEqns and then set y(1)=y_data, but I'm unable to figure it out after several days of fiddling with it.

Akzeptierte Antwort

Star Strider
Star Strider am 22 Sep. 2021
The problem appears to be not allowing the intial conditions to be estimated as parameters as well as the initial value of ‘a_m’.
Also, ‘plot_Times’ must be defined only over the range of times being fitted. If you want to extrapolate to 0 (or any other values), use the deval function.
The only changes I made were to ‘a_m_guess’ in the script, and in the ‘FittingFunction’ function, the ‘InitConds’ assignment and the single parameter ‘a_m(1)’ passed to ‘SysOfDiffEqns’ in the ode45 call.
With only these changes, the code seems to work, and give the appropriate result —
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ]-.35 ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
a_m_guess = [200; 1; 1] ; %init guess for a_m, the parameter to be fitted
a_m = lsqcurvefit(@(am,t)FittingFunction(am,t),a_m_guess,Time_data,Y_data)
Local minimum possible. lsqcurvefit stopped because the size of the current step is less than the value of the step size tolerance.
a_m = 3×1
369.4651 1.0532 4.2399
figure %plot the data as well as the fitted function
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),100) ;
plot(plot_times, FittingFunction(a_m,plot_times),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest')
function F = FittingFunction(a_m,Time_data)
init_slope = 3.39 ; %This is a boundary condition
tspan = Time_data ;
% InitConds = [0 init_slope] ; %y(t=0)=0, y'(t=0)=Us
InitConds = a_m(2:3);
[t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,a_m(1)),tspan,InitConds) ;
F = y_atTdata(:,1) ; %return y only, not y and y'
end
function Y = SysOfDiffEqns(t,y,a_m)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs,
% substitute Y(1) = y
% Y(2) = y'
% makes a system of two first order ODEs:
% Y(1)' = Y(2)
% Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
D = 8.7049 ; % constant input parameter
U_m = 0.95*D ; % costant input parameter
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end
Make appropriate changes to get different results.
.
  8 Kommentare
Christopher Neel
Christopher Neel am 23 Sep. 2021
I figured it out! It finally dawned on me that I can integrate backwards using a regular ode solver by putting in the tspan where it starts at t0 and decreases to zero. That makes all the difference, as it lets me turn my final conditions into initial conditions. In case this helps someone else out in future, here is the full script that does what I wanted. The "extrapolate solution backward" part is the new bit.
clear all
close all
clc
% ydot_o = 3.39 ; %estimate of initial slope
% D_guess = 8.7049 ; %estimate of final slope
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ] ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
%note that the data doesn't start at t=0, so the initial conditions are for *the time at which the data starts*
guess = [1000, 1.5, 3.3, 8.7, 0.94]; % [a_m, y(t_o), y'(t_o), D, U_m_factor] init guess for a_m, y_o, and y'_o, D, Um factor
constrain_parameters = 1 ;
if constrain_parameters == 0
lower_bound = [] ;
upper_bound = [] ;
elseif constrain_parameters == 1
lower_bound = [100, 1, 3, 8, 0.91] ;
upper_bound = [1e4, 2, 4, 9, 0.98] ;
end
return_y_and_ydot = 0 ; % if =0, return only y values (for fitting). If =1, return both y and y' values (for plotting).
options = optimset() ;
[fitted, resnorm, residuals,~,~] = lsqcurvefit(@(fitted,t)FittingFunction(fitted,t,return_y_and_ydot),guess,Time_data,Y_data,lower_bound,upper_bound,options)
%% evaluation and plots
return_y_and_ydot = 1 ;
trans_threshold = 0.99 ;
D = fitted(4) ;
trans_pt = trans_threshold*D ;
figure(1)
subplot(2,1,1) %plot the y(t) data as well as the fitted y(t)
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),1e4) ;
solution = FittingFunction(fitted,plot_times,return_y_and_ydot) ;
plot(plot_times, solution(:,1),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('data and fit vs time'), grid 'on'
subplot(2,1,2) %plot the fitted y'(t) with the transition line
hold 'on'
plot(plot_times, solution(:,2),'-r','DisplayName','fit') ;
line([min(Time_data),max(Time_data)],[trans_pt, trans_pt],'DisplayName','transition')
xlabel('time'), ylabel('Ydot'), legend('Location','NorthWest'), title('ydot vs time'), grid 'on'
figure(2) %plot the residuals
plot(Time_data, residuals,'ko', 'DisplayName','residuals')
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('residuals vs time'), grid 'on'
%% extrapolate solution backwards to t=0
ext_times = linspace(0,min(Time_data),100) ; %extend time to zero
ICs = [fitted(2), fitted(3)] ; %what are the new ICs?
ext_sol = ode23s(@(t,y)SysOfDiffEqns(t,y,fitted),[min(Time_data), 0],ICs) ;
ext_y = deval(ext_sol,ext_times) ;
figure(1)
subplot(2,1,1) % the y(t) plot
plot(ext_times,ext_y(1,:),'-g','DisplayName','extension')
subplot(2,1,2) % the y'(t) plot
plot(ext_times,ext_y(2,:),'-g','DisplayName','extension')
%% functions
function F = FittingFunction(guess,Time_data,return_y_and_ydot)
tspan = Time_data ;
InitConds = [guess(2), guess(3)] ; %y(t=t_o), y'(t=t_o)=init_slope
% [t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
[t_atTdata,y_atTdata] = ode23s(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
if return_y_and_ydot == 0 % used during optimization operation against data
F = y_atTdata(:,1) ; %return y only, not y and y'
elseif return_y_and_ydot == 1 % used to display results
F = y_atTdata ; %return both y and y'
end
end
function Y = SysOfDiffEqns(t,y,guess)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs,
% substitute Y(1) = y
% Y(2) = y'
% makes a system of two first order ODEs:
% Y(1)' = Y(2)
% Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
a_m = guess(1) ;
D = guess(4) ; %let parameter be fitted
U_m_factor = guess(5) ; %let parameter be fitted
% D = 8.7049 ; % constant input parameter
% U_m_factor = 0.94 ; % constant input parameter
U_m = U_m_factor*D ;
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end
Star Strider
Star Strider am 24 Sep. 2021
I dediced to run the code and post it for the benefit of anyone else who encounters this and wants to see the results —
% ydot_o = 3.39 ; %estimate of initial slope
% D_guess = 8.7049 ; %estimate of final slope
Time_data = [ 0.8273 ; 1.1104 ; 1.3940 ; 1.6584 ; 1.8764 ; 2.0173 ; 2.1239 ; 2.2343 ; 2.3498 ] ; % independent variable, experimental data
Y_data = [ 1.5025 ; 2.4650 ; 3.4274 ; 4.3976 ; 5.3600 ; 6.3302 ; 7.2965 ; 8.2590 ; 9.2253 ] ; % dependent variable, experimental data
%note that the data doesn't start at t=0, so the initial conditions are for *the time at which the data starts*
guess = [1000, 1.5, 3.3, 8.7, 0.94]; % [a_m, y(t_o), y'(t_o), D, U_m_factor] init guess for a_m, y_o, and y'_o, D, Um factor
constrain_parameters = 1 ;
if constrain_parameters == 0
lower_bound = [] ;
upper_bound = [] ;
elseif constrain_parameters == 1
lower_bound = [100, 1, 3, 8, 0.91] ;
upper_bound = [1e4, 2, 4, 9, 0.98] ;
end
return_y_and_ydot = 0 ; % if =0, return only y values (for fitting). If =1, return both y and y' values (for plotting).
options = optimset() ;
[fitted, resnorm, residuals,~,~] = lsqcurvefit(@(fitted,t)FittingFunction(fitted,t,return_y_and_ydot),guess,Time_data,Y_data,lower_bound,upper_bound,options)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
fitted = 1×5
1.0e+03 * 1.1093 0.0015 0.0032 0.0087 0.0009
resnorm = 0.0027
residuals = 9×1
0.0131 -0.0226 -0.0019 0.0205 -0.0099 0.0204 -0.0183 -0.0201 0.0187
%% evaluation and plots
return_y_and_ydot = 1 ;
trans_threshold = 0.99 ;
D = fitted(4) ;
trans_pt = trans_threshold*D ;
figure(1)
subplot(2,1,1) %plot the y(t) data as well as the fitted y(t)
hold 'on'
plot(Time_data, Y_data, 'ko', 'DisplayName','data')
plot_times = linspace(min(Time_data),max(Time_data),1e4) ;
solution = FittingFunction(fitted,plot_times,return_y_and_ydot) ;
plot(plot_times, solution(:,1),'-r','DisplayName','fit') ;
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('data and fit vs time'), grid 'on'
subplot(2,1,2) %plot the fitted y'(t) with the transition line
hold 'on'
plot(plot_times, solution(:,2),'-r','DisplayName','fit') ;
line([min(Time_data),max(Time_data)],[trans_pt, trans_pt],'DisplayName','transition')
xlabel('time'), ylabel('Ydot'), legend('Location','NorthWest'), title('ydot vs time'), grid 'on'
figure(2) %plot the residuals
plot(Time_data, residuals,'ko', 'DisplayName','residuals')
xlabel('time'), ylabel('Y'), legend('Location','NorthWest'), title('residuals vs time'), grid 'on'
%% extrapolate solution backwards to t=0
ext_times = linspace(0,min(Time_data),100) ; %extend time to zero
ICs = [fitted(2), fitted(3)] ; %what are the new ICs?
ext_sol = ode23s(@(t,y)SysOfDiffEqns(t,y,fitted),[min(Time_data), 0],ICs) ;
ext_y = deval(ext_sol,ext_times) ;
figure(1)
subplot(2,1,1) % the y(t) plot
plot(ext_times,ext_y(1,:),'-g','DisplayName','extension')
subplot(2,1,2) % the y'(t) plot
plot(ext_times,ext_y(2,:),'-g','DisplayName','extension')
%% functions
function F = FittingFunction(guess,Time_data,return_y_and_ydot)
tspan = Time_data ;
InitConds = [guess(2), guess(3)] ; %y(t=t_o), y'(t=t_o)=init_slope
% [t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
[t_atTdata,y_atTdata] = ode23s(@(t,y)SysOfDiffEqns(t,y,guess),tspan,InitConds) ;
if return_y_and_ydot == 0 % used during optimization operation against data
F = y_atTdata(:,1) ; %return y only, not y and y'
elseif return_y_and_ydot == 1 % used to display results
F = y_atTdata ; %return both y and y'
end
end
function Y = SysOfDiffEqns(t,y,guess)
% y''=a(a_m,y') is the 2nd order diff eq.
% a is a function of a_m() and y'(t)
% need to fit a_m to available pairs of t_data and y_data
% to make system of 1st order deff eqs,
% substitute Y(1) = y
% Y(2) = y'
% makes a system of two first order ODEs:
% Y(1)' = Y(2)
% Y(2)' = a
Y = zeros(2,1) ; % initialize output matrix
a_m = guess(1) ;
D = guess(4) ; %let parameter be fitted
U_m_factor = guess(5) ; %let parameter be fitted
% D = 8.7049 ; % constant input parameter
% U_m_factor = 0.94 ; % constant input parameter
U_m = U_m_factor*D ;
U = y(2) ; % U is dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m)) ; %acceleration function, only valid for U<D
a = max(a,0) ; %force to be at least zero
Y(1) = y(2) ; %1st derivative substitution
Y(2) = a ; %2nd derivative equal to a (which is function of y')
end
I congratulate on extending the original idea!
.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Contour Plots finden Sie in Help Center und File Exchange

Produkte


Version

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by