Solving second order ODEs with ANN

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

[Murat Balc?]

Hi Dimitris, aren't the following pieces of code identical to each other?

% 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"));

Answer 1

[Ranjeet]

Hi Dimitris,

I assume that the code you provided has been adopted from the following resource Solve Ordinary Differential Equation Using Neural Network.

The resource guides on training a neural network to estimate solution of a first order differential equation. The important point to be noted is the trained network estimates the analytical solution of the taken example ().

For the case of second order differential equation, we get the analytical solution in many forms including , etc. Here, are found from characteristic equation and are found from the initial conditions (Analytical solution exist in other forms as well).

Changes suggested to adopt the example for second order differential equation are:

1. Modify the modelGradients function to adopt the loss calculation based on second order differential equation taken to solve.
2. Test by increasing number of hidden layers as the analytical solution is now more complex.
3. Try changing other hyperparameters based on the obtained test results.

[Dimitris Kastoris]

thank you very much!