Can i use this method for second order ODEs and is this modification on modelGradients correct?
x = linspace(0,1,10000)';
inputSize = 1;
layers = [
featureInputLayer(inputSize,Normalization="none")
fullyConnectedLayer(10)
sigmoidLayer
fullyConnectedLayer(1)
sigmoidLayer];
dlnet = dlnetwork(layers);
numEpochs = 15;
miniBatchSize =100;
initialLearnRate = 0.1;
learnRateDropFactor = 0.3;
learnRateDropPeriod =5 ;
momentum = 0.9;
icCoeff = 7;
ads = arrayDatastore(x,IterationDimension=1);
mbq = minibatchqueue(ads,MiniBatchSize=miniBatchSize,MiniBatchFormat="BC");
figure
set(gca,YScale="log")
lineLossTrain = animatedline(Color=[0.85 0.325 0.098]);
ylim([0 inf])
xlabel("Iteration")
ylabel("Loss (log scale)")
grid on
velocity = [];
iteration = 0;
learnRate = initialLearnRate;
start = tic;
% Loop over epochs.
for epoch = 1:numEpochs
% Shuffle data.
mbq.shuffle
% Loop over mini-batches.
while hasdata(mbq)
iteration = iteration + 1;
% Read mini-batch of data.
dlX = next(mbq);
% Evaluate the model gradients and loss using dlfeval and the modelGradients function.
[gradients,loss] = dlfeval(@modelGradients3, dlnet, dlX, icCoeff);
% Update network parameters using the SGDM optimizer.
[dlnet,velocity] = sgdmupdate(dlnet,gradients,velocity,learnRate,momentum);
% To plot, convert the loss to double.
loss = double(gather(extractdata(loss)));
% Display the training progress.
D = duration(0,0,toc(start),Format="mm:ss.SS");
addpoints(lineLossTrain,iteration,loss)
title("Epoch: " + epoch + " of " + numEpochs + ", Elapsed: " + string(D))
drawnow
end
% Reduce the learning rate.
if mod(epoch,learnRateDropPeriod)==0
learnRate = learnRate*learnRateDropFactor;
end
end
ModelGradients
function [gradients,loss] = modelGradients2(dlnet, dlX, icCoeff)
y = forward(dlnet,dlX);
% Evaluate the gradient of y with respect to x.
% Since another derivative will be taken, set EnableHigherDerivatives to true.
dy = dlgradient(sum(y,"all"),dlX,EnableHigherDerivatives=true);
% Define ODE loss.
eq = dy + y/5 - exp(-(dlX / 5)) .* cos(dlX);
% Define initial condition loss.
ic = forward(dlnet,dlarray(0,"CB")) - 0 ;
% Specify the loss as a weighted sum of the ODE loss and the initial condition loss.
loss = mean(eq.^2,"all") + icCoeff * ic.^2;
% Evaluate model gradients.
gradients = dlgradient(loss, dlnet.Learnables);
end
New modelGradients
function [gradients,loss] = modelGradients4(dlnet, dlX, icCoeff)
y = forward(dlnet, dlX);
% Evaluate the gradient of y with respect to x.
dy = dlgradient(sum(y, "all"), dlX, 'EnableHigherDerivatives', true);
d2y = dlgradient(sum(dy, "all"), dlX, 'EnableHigherDerivatives', true);
% Define ODE loss.
eq = d2y + 1/5*dy + y + 1/5*exp(-dlX/5).*cos(dlX);
% Compute the initial conditions directly within this function
yAtZero = forward(dlnet, dlarray(0,"CB"));
% Use the gradient dy and evaluate it at x = 0
dyAtZero = forward(dlnet, dlarray(0,"CB"));
% Initial condition errors
icErrorY = yAtZero ; % because y(0) = 0
icErrorDY = dyAtZero - 1; % because y'(0) = 1
% Combine the ODE loss and the initial condition errors.
loss = mean(eq.^2, "all") + icCoeff * (icErrorY.^2 + icErrorDY.^2);
% Evaluate model gradients.
gradients = dlgradient(loss, dlnet.Learnables);
end
