Hello,
I am trying to find the optimal parameters for my first order ODE to fit a dataset.
My ODE takes the following form: ds/dt= (alfa * temp – s(t)) * (1/tau)
The parameters that I need to find are "alfa" and "tau" whereas "temp" (temperature) depends on time but I have no equation to represent it, just a dataset of its evolution over time during the same timeframe of "ds/ dt" evolution.
I would appreciate any feedback on how to get "temp" into the code as my working version is not working. Please see below the code and its errors.
Many thanks in advance for any feedback.
CODE:
%Initial data
load rcp85_expansionmidmatlab.txt;load rcp85_temperaturemidmatlab.txt;time= rcp85_temperaturemidmatlab(:,1)-2006+1;temp= rcp85_temperaturemidmatlab(:,2);sealevel=rcp85_expansionmidmatlab(:,2);plot(time,sealevel)%ODE information
tSpan = [1 95];z0 = [0.0118];tempi=temp(tSpan);%Initial guess
alfa=0.2;%lower range from Mengel paper
tau=82;%lower range from Mengel paperODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau,tempi),tSpan,z0); % Run the ODE
simsealevel = deval(ODE_Sol, time); % Evaluate the solution at the experimental time steps
hold onplot(time, simsealevel, '-r')%% Set up optimization
myObjective = @(x) objFcn(x,time,sealevel,tSpan,z0);lb = [0.2,82];ub = [0.63,1290];bestalfatau = lsqnonlin(myObjective, alfa,tau, lb, ub);%% Plot best result
ODE_Sol = ode45(@(t,z)updateStates(t,z,bestalfatau,tempi),tSpan, z0);bestsealevel = deval(ODE_Sol, time);plot(time, bestsealevel, '-g')legend('IPCC Data','Initial Param','Best Param');function f = updateStates(t,z,alfa,tau,tempi)f = (alfa.*tempi-z)*(1/tau);function cost= objFcn (x,time,sealevel,tSpan,z0)ODE_Sol = ode45(@(t,z)updateStates(t,z,x),tempi, tSpan, z0);simsealevel = deval(ODE_Sol, time);cost = simsealevel-sealevel;
ERRORS:
Error using odearguments (line 95)@(T,Z)UPDATESTATES(T,Z,ALFA,TAU,TEMPI) returns a vector of length 2, but the length of initial conditions vector is 1.The vector returned by @(T,Z)UPDATESTATES(T,Z,ALFA,TAU,TEMPI) and the initial conditions vector must have the same numberof elements.Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);Error in odeparam1 (line 21)ODE_Sol= ode45(@(t,z)updateStates(t,z,alfa,tau,tempi),tSpan,z0); % Run the ODE
Best Answer