HW5 Ex1

import numpy as np
import torch
import torch.nn as nn
from torch import optim
import matplotlib.pyplot as plt

# υπολογισμός παραγώγου του  u
def deriv_u(u, x):
    from torch import autograd
    grad, = autograd.grad(outputs = u, inputs = x, grad_outputs = torch.ones_like(u), create_graph = True)
    return grad_u

# υπολογισμός παραγώγου του  w
def deriv_w(w, x):
    from torch import autograd
    grad, = autograd.grad(outputs = w,inputs = x, grad_outputs = torch.ones_like(w) , create_graph = True)
    return grad_w

neurons = 50
netu = nn.Sequential(nn.Linear(1, neurons), nn.Sigmoid(), nn.Linear(neurons,1,bias = False))

neurons = 50
netw = nn.Sequential(nn.Linear(1, neurons), nn.Sigmoid(), nn.Linear(neurons,1,bias = False))

# πεδίο ορισμού του t
x = torch.linspace(0.0, 2*np.pi, 100, requires_grad = True).reshape(-1,1)

loss_fun = loss_u + loss_w

optimizer = optim.Adam(list[netu.parameters()+list(netw.parameters())], lr = 0.01)

# u

epochs = 1000
loss_hist_u = []
u_hist = []
for i in range (epochs):
    optimizer.zero_grad()
    u = x + x*netu(x)
    dudx = deriv_u(u,x)
    loss_u = nn.MSELoss()
    loss_w = nn.MSELoss()
    loss = loss_fun(dudx, w(x))
    loss_hist_u.append(loss.item())
    loss.backward()
    optimizer.step()
    u_hist.append(u.detach())
    
# εύρεση της λύσης του  u  με κλασική μέθοδο

from scipy.integrate import solve_ivp
t = np.linspace(0, 2*np.pi, 100)
def ode1(x,u):
    return w(x)

solu = solve_ivp(ode1, [0, 2*np.pi], [1], t_eval=t, method = 'RK45')

# w

epochs = 1000
loss_hist_w = []
w_hist = []
for i in range (epochs):
    optimizer.zero_grad()
    w = (1-x)+x*netw(x)
    dwdx = deriv_w(w,x)
    loss_u = nn.MSELoss()
    loss_w = nn.MSELoss()
    loss = loss_fun(dwdx, -u(x))
    loss_hist_w.append(loss.item())
    loss.backward()
    optimizer.step()
    w_hist.append(w.detach())
    
# εύρεση της λύσης του  w  με κλασική μέθοδο

from scipy.integrate import solve_ivp
t = np.linspace(0, 2*np.pi, 100)
def ode2(x,w):
    return -u(x)

solw = solve_ivp(ode2, [0, 2*np.pi], [1], t_eval=t, method = 'RK45')
    
# plot

plt.figure()
plt.plot(x.detach(), u.detach(), 'b', label = "neural network")
plt.plot(t,solu.y[0], 'r--', label = "numerical solution with RK 45")
plt.legend()
plt.show()

plt.figure()
plt.plot(loss_hist_u, 'b')
plt.show()

plt.figure()
plt.plot(x.detach(), w.detach(), 'b', label = "neural network")
plt.plot(t,solw.y[0], 'r--', label = "numerical solution with RK 45")
plt.legend()
plt.show()

plt.figure()
plt.plot(loss_hist_w, 'b')
plt.show()

# HOMEWORK 5 ΑΣΚΗΣΗ 1