1
$\begingroup$

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;
$\endgroup$
  • $\begingroup$ What happens if you set InterpolationOrder to 1: xInt = Interpolation[XTbl, IterpolationOrder->1]? $\endgroup$ – Quantum_Oli Oct 26 '16 at 16:20
  • $\begingroup$ @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 $\endgroup$ – Serge A. Sergeev Oct 26 '16 at 17:59
  • $\begingroup$ @andre Don't quite understand what you mean. $\endgroup$ – Serge A. Sergeev Oct 26 '16 at 18:01
  • $\begingroup$ @andre I've found it - while copying here made syntax error. But the question remains. $\endgroup$ – Serge A. Sergeev Oct 26 '16 at 18:25
  • $\begingroup$ @andre I added initial conditions for the system at the end of the post. $\endgroup$ – Serge A. Sergeev Oct 26 '16 at 18:49
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.