NDSolve for non-analytical forms

Question

I came across a problem of coupled differential equations of a non-analytical results. It follows that NDsolve requires an evaluated form of equations to be fed into it for a differential equation to process.But, my functions can't be evaluated analytically. I want to know if there is any possible way to tackle the problem. I have created a sample problem here:

First I create a Complex system just to confirm I will get non-analytical solutions:

Mat[n_, x_, y_] := SparseArray[{Band[{1, 1}, {n, n}] -> {x^3, x I + 5 y^2 + 4, Sqrt[x]},
                   Band[{1, 2}, {n, n}] -> {y^3, Sqrt[x + I x^3 - y^2 + 4], 
                   Sqrt[x - y^2]}, Band[{3, 1}, {n, n}] -> {I x^3, x + y x^2 + 4, Sqrt[x + y^2]}}]
eval[n_, x_, y_] := Eigensystem[Mat[n, x, y]][[1]]
evec[n_, x_, y_] := Eigensystem[Mat[n, x, y]][[2]]

Secondly, I will form some time-dependent functions so as to cover the complexity of my problem:

val = D[Mat[8, x, y], x];
x2[x1_, t_] := x1 + t^2
y2[y1_, t_] := y1 + t
T1[x1_, y1_] := {I Conjugate[#], #} &@(Conjugate[evec[8, x1, y1][[6]]].SparseArray[
                ArrayRules[val] /. {x -> x1, y -> y1}, Dimensions[val]].evec[8, x1, y1][[7]] // N)
T2[t_] := {t, t^2}
T3[x1_, y1_, t_] := T1[x2[x1, t], y2[y1, t]]
T4[x1_, y1_, t_] := Re[eval[8, x2[x1, t], y2[x1, t]][[7]]] - Re[eval[8, x2[x1, t], y2[x1, t]][[8]]]//N

Last, I will list them in the required form of equation and initial conditons for processing:

t0=-5;
eqns[x1_, y1_, t_] := {A1'[t] == (T2[t].T3[x1, y1, t]) A2[t], A2'[t] == A1[t] (T2[t].Conjugate[T3[x1, y1, t]]),
                      i'[t] == T4[x1, y1, t],A1[t0] == 0, A2[t0] == 1, i[t0] == 0}
sol1 = ParametricNDSolve[eqns[x1, y1, t], {A1[t], A2[t]}, {t, t0, 5}, {x1, y1}]

As you see, say, eqns[x1, y1, t] can't be evaluated unless you provide numerical values of all the parameters. How do we solve the equations in that case. I would be grateful for your help.

(Note: this is a sample just to reflect my problem, feel free to make reasonable changes)

eqns also has a function i. Have you defined it or should it be solved for? — Natas, Jun 19 '20 at 11:05
It’s a way of writing if there is an integral in your calculations . I want to evaluate integral of ‘T4’ alongside. So, Yes, ‘i’ should be solved for. I have provided initial conditions for it. — Rupesh, Jun 19 '20 at 11:17
I am really not sure what you are trying to do here. Can you simplify your example to something much simpler? Are you looking for the NumericQ functionality, as described here? — MarcoB, Jun 19 '20 at 15:25
Sure, Let's just solve for a parameter containing T4. Much more simplified, how can i solve this: ParametricNDSolve[ i'[t] == T4[x1, y1, t], i[t0] == 0, {i[t]}, {t, t0, 5}, {x1, y1}]. — Rupesh, Jun 19 '20 at 17:18
@MarcoB All I am trying is to solve a differential equation using NDsolve. The equations can only be numerically determined. Above, In the comment I have presented simplified version of it. Please do mention if you want any further info — Rupesh, Jun 19 '20 at 17:55

I.M. · Accepted Answer · 2020-06-21T09:11:22.550

5

Your equations are of the form x' = f(x), x = {x1,...,xn}, flow part f(x) can be defined implicitly, i.e. as a black box outside of NDSolve.

This explicit form:

sol = NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2x[t]-y[t]^3,x[0]==y[0]==1},{x,y},{t,20}]
ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,20}]

can be converted to:

ClearAll[flow] ;
flow[x_,y_] := flow[x,y] = {-y-x^2,2 x-y^3} ;
ClearAll[fx,fy] ;
fx[arg__?NumericQ] := flow[arg][[1]] ;
fy[arg__?NumericQ] := flow[arg][[2]] ;
sol = NDSolve[{x'[t] == fx[x[t],y[t]] ,y'[t]==fy[x[t],y[t]],x[0]==y[0]==1},{x,y},{t,20}] ;
ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,20}]

Edit

In your case:

ClearAll[flow] ;
flow[x1_, y1_, t_,A1_,A2_,i_] :=  flow[x1,y1,t,A1,A2,i] = {
    (T2[t].T3[x1, y1, t]) A2,
    A1 (T2[t].Conjugate[T3[x1, y1, t]]),
    T4[x1, y1, t]
} ;
ClearAll[f1,f2,f3] ;
f1[arg__?NumericQ] := flow[arg][[1]] ;
f2[arg__?NumericQ] := flow[arg][[2]] ;
f3[arg__?NumericQ] := flow[arg][[3]] ;
t0=-5;
(* return only i *)
sol = ParametricNDSolveValue[
    {
        A1'[t] == f1[x1,y1,t,A1[t],A2[t],i[t]],
        A2'[t] == f2[x1,y1,t,A1[t],A2[t],i[t]],
        i'[t] == f3[x1,y1,t,A1[t],A2[t],i[t]],
        A1[t0] == 0, A2[t0] == 1, i[t0] == 0
    },
    i,
    {t, t0, 5},
    {x1, y1}
] 
sol[1,1]

edited Jun 21 '20 at 09:11

answered Jun 21 '20 at 06:03

I.M.

2,926
1
13
18

Thank you. When you defined flow[x,y] = {-y-x^2,2 x-y^3} that still becomes an analytical representation. Whereas, if you see in mine it has a set delayed. Could you just possibly obtain ParametricNDSolve[ i'[t] == T4[x1, y1, t], i[t0] == 0, {i[t]}, {t, t0, 5}, {x1, y1}] from my example. – Rupesh Jun 21 '20 at 06:21
@Rupesh, you can define flow[x_,y_] := flow[x,y] = whatever, flow[x_,y_] := flow[x,y] part is a kind of memorization, since both fx and fy evaluate flow with the same argument – I.M. Jun 21 '20 at 06:33
I understand what you are saying but since there are many intermediary steps from defining a flow to NDSolve in my case, how do you go on about it. Plus, I am trying to find a time evolution of a separate function but not x an y itself. I am afraid that I can't make a resemblance between your suggestion and my problem. I would be glad if you could take my own example and at least find an evolution for i[t]. You can change the time intervals or I expect the solution parametric in terms of x and y – Rupesh Jun 21 '20 at 06:54
@Rupesh, see the edit. – I.M. Jun 21 '20 at 09:12
@ I.M I want to thank you a lot and appreciate your time and patience. I was totally unaware of this approach and now it does make sense to me. I tried fitting your first example to my case and messed up, it was bad on my part. Sorry if I behaved impolitely. – Rupesh Jun 21 '20 at 12:38
@Rupesh, np, glad to help – I.M. Jun 21 '20 at 13:20
Hi I.M I am reaching you out because my actual problem still isn't solved. I tried all the available methods but the NDSolve just gets stuck. I have tried the regular approach and the one you suggested. I would be really grateful if you could look into my problem when you have time. Here is the link to the notebook – Rupesh Jul 06 '20 at 00:25
@Rupesh, rewriting of equations seems to work – I.M. Jul 06 '20 at 03:06
Thank you for your help. On the side note the indexing way which I wrote earlier works too, maybe Through is the key.Here is a pic of using A1[1,t] and A11.Another thing, the runtime is almost 10 mins which I feel is very huge considering the computation involved, did you have a similar experience?. Thank you again for your help – Rupesh Jul 06 '20 at 05:19
@Rupesh, indeed, this indexing does work, still, it is confusing. Might be so for NDSolve[] too, because your eqs look like pdes. Tested in wolframclould, took less than 1m to compute for one fixed value of qx. I've changed indexing and added memorization sorted[n_,kx_]:=sorted[n,kx] = .... Through should not effect performance here. – I.M. Jul 06 '20 at 05:39
Thanks again. It's an ODE in time which does contain some parametric partial derivatives. I agree for a two-variable scenario indexing isn't required but If we set Range to different, then you will have multiple variables and different sets of initial conditions. see fig. For 1 set of initial condition initials[[1]], our ODE will look as pic. For this, it's stuck. Would you feel free to try it for a higher range scenario in your cloud computation? .Modified file. – Rupesh Jul 06 '20 at 13:53
@Rupesh, it takes about 2s per evaluation with basic wolframcloud plan, you can try to adjust the integration method or optimize the computation of functions involved, file. – I.M. Jul 06 '20 at 19:02
You are amazing and really good at scripting. I will go through it and hopefully it will be solved. One final request, would you remove the dropbox link as it's my project file. Also, I will cite you if I get anything publishable, I'll let you know. Thank you so much for your help. – Rupesh Jul 06 '20 at 20:35

NDSolve for non-analytical forms

1 Answers1

Edit