# How to interpolate data on unstructured grid?

I have the Hamilton system in 3D case

    HamSysFull[Psi_, Theta_,  alpha_] :=
NDSolve[{X'[t] ==
Hp[X[t], Y[t],  Z[t], P[t], Q[t], Omega[t]],
Y'[t] == Hq[X[t], Y[t],  Z[t], P[t], Q[t], Omega[t]],
Z'[t] ==
HOmega[X[t], Y[t],  Z[t], P[t], Q[t],  Omega[t]],
P'[t] == -Hx[X[t], Y[t],  Z[t], P[t], Q[t], Omega[t]],
Q'[t] == -Hy[X[t], Y[t],  Z[t], P[t],
Q[t],  Omega[t]], Omega'[t] == -Hz[X[t], Y[t], Z[t], P[t], Q[t],   Omega[t]],
X[0] == Xst[alpha, Psi, Theta],
Y[0] == Yst[alpha, Psi, Theta],
Z[0] == Zst[alpha, Psi, Theta],
P[0] == Pst[Psi, Theta], Q[0] == Qst[Psi, Theta],
Omega[0] == Omegast[Psi, Theta]},  {X,
Y, Z, P, Q, Omega}, {t, 0, NumT},
StepMonitor :> { Sow[{{t, Psi, Theta}, X[t]}, x],
Sow[{{t, Psi, Theta}, Y[t]}, y],
Sow[{{t, Psi, Theta}, Z[t]}, z],
Sow[{{t, Psi, Theta}, P[t]}, p],
Sow[{{t, Psi, Theta}, Q[t]}, q],
Sow[{{t, Psi, Theta}, Omega[t]}, Omega]}]


After I solve if for different angles psi and theta

angleBound = Pi/4;
htheta = 2 angleBound/10;
hpsi = Pi/10;
fullSol = Table[Reap[HamSysFull[Psi, Theta, 0];], {Psi, 0,
2 Pi, hpsi}, {Theta, -angleBound, angleBound, htheta}];


And for each solution I build the list for interpolation

 XTbl = Flatten[Table[fullSol[[i, j, 2, 1]],
{i, 1, Length[fullSol]}, {j,1,Length[fullSol[[i]]]}], 2];


But the command

xInt = Interpolation[XTbl]


returns the error

 Interpolation::udeg: Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will be reduced to 1. >>


My question is how to interpolate such list of data. As I assumed it is because the unstructured grid by time. Making the structured grid is not a good option because it takes a lot of evalution.

Initial conditions for the system

Hdepth = 4000;
Epsilon = 0.00737;
A = 1500;
beta = 2;
z0 = 1300;
Zf[z_] := beta (z - z0)/z0;
c[z_] := A (1 + Epsilon (Zf[z] - 1 + Exp[-Zf[z]]));

H[X_, Y_,  Z_, P_, Q_,  Omega_] :=  c[-Z] Sqrt[P^2 + Q^2 + Omega^2]
Hp[X_, Y_,  Z_, P_, Q_,  Omega_] := D[H[X, Y, Z, P, Q, Omega], P];
Hq[X_, Y_,  Z_, P_, Q_,  Omega_] := D[H[X, Y, Z, P, Q, Omega], Q];
HOmega[X_, Y_,  Z_, P_, Q_, Omega_] := D[H[X, Y, Z, P, Q, Omega], Omega];
Hx[X_, Y_,  Z_, P_, Q_, Omega_] := D[H[X, Y, Z, P, Q, Omega], X];
Hy[X_, Y_,  Z_, P_, Q_, Omega_] := D[H[X, Y, Z, P, Q, Omega], Y];
Hz[X_, Y_,  Z_, P_, Q_, Omega_] := D[H[X, Y, Z, P, Q, Omega], Z];

Clear[x, z, p, Omega, Xpsi, Zpsi, Ppsi, Omegapsi, solution, t];

stDepth = -2000;  (* depth by z *)
stShiftx = 0;  (* shifting by x *)
stShifty = 0;  (* shifting by  *)
stPoint = {stShiftx, stShifty, stDepth}; (* starting conditions *)
alpha = 0; (* coordinate on manifold *)

Pst[Psi_, Theta_] := Cos[Psi] Cos[Theta];
Qst[Psi_, Theta_] := Sin[Psi] Cos[Theta];
Omegast[Psi_, Theta_] := Sin[Theta];
Xst[Alpha_ , Psi_, Theta_] := stPoint[[1]];
Yst[Alpha_ , Psi_, Theta_] := stPoint[[2]];
Zst[Alpha_ , Psi_, Theta_] := stPoint[[3]];

NumT = 10;

• What happens if you set InterpolationOrder to 1: xInt = Interpolation[XTbl, IterpolationOrder->1]? Oct 26, 2016 at 16:20
• @Quantum_Oli It gives another series of errors about the grid: Interpolation::femimq: The element mesh has insufficient quality of -5.22882*10^-16. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements. >>8 Oct 26, 2016 at 17:59
• @andre Don't quite understand what you mean. Oct 26, 2016 at 18:01
• @andre I've found it - while copying here made syntax error. But the question remains. Oct 26, 2016 at 18:25
• @andre I added initial conditions for the system at the end of the post. Oct 26, 2016 at 18:49

## 1 Answer

this requires no more computation, remove the StepMonitor and move the Sow to a table evaluated after NDSolve has finished:

HamSysFull[Psi_, Theta_, alpha_] :=
Module[{res},
res = First@
NDSolve[{X'[t] == Hp[X[t], Y[t], Z[t], P[t], Q[t], Omega[t]],
Y'[t] == Hq[X[t], Y[t], Z[t], P[t], Q[t], Omega[t]],
Z'[t] == HOmega[X[t], Y[t], Z[t], P[t], Q[t], Omega[t]],
P'[t] == -Hx[X[t], Y[t], Z[t], P[t], Q[t], Omega[t]],
Q'[t] == -Hy[X[t], Y[t], Z[t], P[t], Q[t], Omega[t]],
Omega'[t] == -Hz[X[t], Y[t], Z[t], P[t], Q[t], Omega[t]],
X[0] == Xst[alpha, Psi, Theta], Y[0] == Yst[alpha, Psi, Theta],
Z[0] == Zst[alpha, Psi, Theta], P[0] == Pst[Psi, Theta],
Q[0] == Qst[Psi, Theta], Omega[0] == Omegast[Psi, Theta]}, {X,
Y, Z, P, Q, Omega}, {t, 0, NumT}];
Table[Sow[{{t, Psi, Theta}, (X[t] /. res)}, x], {t, 0, NumT, NumT/100}];
Table[Sow[{{t, Psi, Theta}, (Y[t] /. res)}, y], {t, 0, NumT, NumT/100}]]


now with the rest of your code unchanged, Interpolation[XTbl] works fine.