1
$\begingroup$

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)

$\endgroup$
  • 1
    $\begingroup$ eqns also has a function i. Have you defined it or should it be solved for? $\endgroup$ – Natas Jun 19 '20 at 11:05
  • $\begingroup$ 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. $\endgroup$ – Rupesh Jun 19 '20 at 11:17
  • 1
    $\begingroup$ 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? $\endgroup$ – MarcoB Jun 19 '20 at 15:25
  • $\begingroup$ 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}]. $\endgroup$ – Rupesh Jun 19 '20 at 17:18
  • $\begingroup$ @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 $\endgroup$ – Rupesh Jun 19 '20 at 17:55
5
+50
$\begingroup$

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]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Rupesh Jun 21 '20 at 6:21
  • $\begingroup$ @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 $\endgroup$ – I.M. Jun 21 '20 at 6:33
  • $\begingroup$ 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 $\endgroup$ – Rupesh Jun 21 '20 at 6:54
  • $\begingroup$ @Rupesh, see the edit. $\endgroup$ – I.M. Jun 21 '20 at 9:12
  • $\begingroup$ @ 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. $\endgroup$ – Rupesh Jun 21 '20 at 12:38

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.