Pull Request: Implement PINN Solver for Timoshenko Beam Theory with Two-Stage Optimization

[DenizCanELCI]

Summary
This PR replaces the single 4th-order Euler-Bernoulli beam formulation with a coupled system of two 2nd-order differential equations based on Timoshenko beam theory. It captures shear deformation kinematics alongside bending deflection and addresses gradient stalling by introducing an Adam-to-L-BFGS training pipeline.

Motivation & Physical Context
In thick beam regimes (deep beams where h/L is non-negligible), classic Euler-Bernoulli theory underestimates total compliance by ignoring shear deformation.

Moving to Timoshenko theory allows independent tracking of vertical deflection (w) and cross-section rotation (phi).

Formulating the problem as two coupled 2nd-order PDEs reduces the autodiff computational graph depth compared to 4th-order derivatives, mitigating vanishing gradient artifacts during backpropagation.

Key Changes
Governing Equations & Multi-Output Architecture:

Updated the neural network output dimension to 2: y_pred = [w, phi].

Implemented coupled equilibrium residuals via dde.grad.jacobian and dde.grad.hessian:

Shear equilibrium: kGA * (w'' - phi') + q = 0

Moment equilibrium: EI * phi'' + kG*A * (w' - phi) = 0

Boundary Conditions (Cantilever Setup):

Clamped Root (x = 0): Enforced Dirichlet BCs on both components (w(0) = 0, phi(0) = 0). Total centerline slope w'(0) = V(0) / (kGA) != 0 due to root shear reaction.

Free Tip (x = 1): Enforced natural boundary conditions via OperatorBC with explicit tensor indexing:

Moment: M(1) = 0 --> d(phi)/dx = 0

Shear Force: V(1) = 0 --> dw/dx - phi = 0

Optimization Strategy:

Resolved the training plateau where the first-order Adam optimizer stalled on natural boundary condition gradients (L2 error plateaued around ~70%).

Added a two-phase optimization schedule: initial rough convergence using Adam followed by quasi-Newton L-BFGS fine-tuning.

Analytical Ground Truth & Verification:

Corrected the shear contribution sign in the closed-form benchmark function:

Python
w_total = w_bending + (1.0 / GA_s) * (0.5 * x**2 - x)

[echen5503]

keep in mind that the repository main language is english (secondary chinese), and please also delete the loss files.

[DenizCanELCI]

Thank you @echen5503 . I updated some comments in my original language. I also deleted the .dat files. Best regards.