ode45 array issues

1 visualización (últimos 30 días)
Jorge Cossio
Jorge Cossio el 16 de Abr. de 2019
Comentada: Torsten el 16 de Abr. de 2019
Running into issues with the familiar "Index exceedes array bounds" while using ode45.
Trying to figure out the value for k in the above code based on a system below and the dataset 'Data.MAT' consisting of 24 concentration and time measurements.
Can't figure out what's going on, I've tried defining my equation in a separate function (gave up on that) and instead followed the example on this Mathworks ode 45 example. I've used the same method before but can't understand why it's not working this time.
clc
clear all
close all
%% Loading and viewing data
Data = load('ProbData_hw3.MAT');
tspan = Data.time;
Ca_data = Data.Soln;
plot(tspan,Ca_data,'ko ')
xlabel('time')
ylabel('Measurement')
%% Setting up
y0 = Ca_data(1);
learn_rate = .001;
k_started = [1 1.001];
%get things going to get dJ/dk
% Get J1
k = k_started(1);
b=sqrt(1.2*9.81*k);
syms y(t)
[V] = odeToVectorField(diff(y,2) == (-k/1)*(9.81/k +b*y) - diff(y));
f = matlabFunction(V,'vars', {'t','Y'});
[t, y] = ode45(f, [0 5], y0);
y_hat = interp1(t,y,tspan);
for n = 1:24
tmp(n) = (y_hat(n)-Ca_data(n))^2;
end
J1 = 1/n*sum(tmp);
for i = 1:500
if i == 1
k = k_started(2);
Jold = J1;
k_old = k_started(1);
else
k = k_new;
Jold = Jnew;
end
[t, y] = ode45('ConcODE', [0 5], y0, []);
y_hat = interp1(t,y,tspan);
for n = 1:24
tmp(n) = (y_hat(n)-Ca_data(n))^2;
end
Jnew = 1/n*sum(tmp);
dP = -learn_rate*(Jnew-Jold)/(k-k_old); % k(+1) = k - LR(dJ/dK)
k_old = k;
k_new = k+dP;
Ki(i) = k_new;
Ji(i) = Jnew;
if dP < 1e-6
break
end
end
hold on
plot(t,y)
figure, plot(Ki,Ji)
  2 comentarios
Walter Roberson
Walter Roberson el 16 de Abr. de 2019
Unfortunately we cannot test without the .mat
Torsten
Torsten el 16 de Abr. de 2019
  1. y0 must be a 2-element vector
  2. The ODE you gave in the text and the one used in the code seem to differ.
Try
to solve your ODE analytically first and then use lsqcurvefit to fit k to your data.

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by