Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica: NDSolve with InterpolatingFunction as initial conditions

  1. Oct 26, 2011 #1
    The questions arises because I want to use the solution of one PDE as initial condition to solve another. Then using NDSolve, the first solution is given by InterpolatingFunction. I tried, it takes the whole afternoon, almost kills my computer but still not giving any result. Another computer simply refused to work after sometime, telling me memory is not enough. My question is:

    1. Can InterpolatingFunction be used as an initial condition in NDSolve?

    2. Is there any way to bypass InterpolatingFunction? For example, make a list of values out of InterpolatingFunction, and then solve PDE based on this list of initial conditions, just like Plot and ListPlot?
  2. jcsd
  3. Oct 26, 2011 #2
    Can you post some example code? Not necessarily the whole thing you're working on, but a simplified example that contains the essence of the problem...
  4. Oct 26, 2011 #3
    Here is an example interpolating function and a table created from it.


    Out[1]= {1/2, 1/3, 1/4, 1/5, 2/15}
  5. Oct 26, 2011 #4
    So long as what you plug in for your initial conditions is a number you will be fine, i.e. pre-evaluate the interpolating function before you try ndsolve (use a loop or table or something).
  6. Oct 27, 2011 #5
    Hi, All, here is my example:
    c = 29.979;
    T = 500;
    tm = -250;
    σ = 100;

    solS=NDSolve[{∂tsS[z,t]==-(γ+I β*z)sS[z,t]-I*w*aS[z,t], (∂zaS[z,t]+1/c ∂t aS[z,t])==-I*h*k sS[z,t], sS[z,-T]==0, aS[z,-T]==E^((-T-tm)2/2σ2), aS[L0,t]==E^((t-tm)2/2σ2),{sS,aS},{z,L0,L1},{t,-T,0}];


    solS=NDSolve[{∂tsR[z,t]==-(γ-I β z)sR[z,t]-I*w*aR[z,t], (∂zaR[z,t]+1/c ∂t aR[z,t])==-I*h*k sR[z,t], sR[z,0]==iniR[z,0], aR[z,0]==0, aR[L0,t]==0,{sR,aR},{z,L0,L1},{t,0,T}]
    Last edited: Oct 27, 2011
  7. Oct 27, 2011 #6
    FunkyDwarf: The initial condition is not just A number, it should be a function. But the analytical form of the function is unknown, because it is the solution of another PDE. See my example above.
  8. Oct 27, 2011 #7

    Bill: Yes, one can make a list like this. But the question is how to make it into an initial condition of NDSolve?
  9. Oct 27, 2011 #8
    Simon: Example is given in 5th floor.
  10. Oct 27, 2011 #9
    You need to first get just the first one working and that code above is not usable without a lot of alterations because of the special characters in it. To post Mathematica code in here, first select it in the notebook, then do a Cell/Convert To/Raw input form. Then cut and paste the raw input form it here encolsed in code tags selected from the menu above.

    Also, your first one needs six conditions, three for each function. You need the initial condition say sS(z,0)=myf(z), then the two boundary conditions like sS(-l0,t)=0 and ss(l1,t)=something else. Do that for each function.

    Don't start user-defined variables with upper case letter. Can conflict with built-in functions.

    Ok, then just get the first one working where Mathematica reports Interpolation functions for both the variables sS and aS, then we can work from there.
  11. Oct 27, 2011 #10
    jackmell: Thanks for telling me. Well, I tried to paste but it looked ugly. I'll try more. And the first PDE indeed works, though it has only three boundary/initial conditions. I'll show later.
  12. Oct 27, 2011 #11
    Working Mathematica Note:

    Thanks for telling me how to insert code. Typing code directly was killing me.

    Code (Text):
    \[Gamma] = 10^(-5);
    T = 500;
    tm = -250;
    \[Sigma] = 100;
    L = 1;
    L0 = -(L/2);
    L1 = L0 + L;
    c = 29.979;
    \[Omega] = 1.32949/10^8;
    \[Eta] = 1539.1;
    \[Mu] = 6.27;
    \[Beta] = 0.1334;
    k = 40895.3;
    solS = NDSolve[{D[sS[z, t], t] == (-(\[Gamma] + I*\[Beta]*z))*sS[z, t] - I*\[Omega]*aS[z, t],
        D[aS[z, t], z] + (1/c)*D[aS[z, t], t] == (-I)*\[Eta]*k*sS[z, t], sS[z, -T] == 0,
        aS[z, -T] == E^(-((-T - tm)^2/(2*\[Sigma]^2))),
        aS[L0, t] == E^(-((t - tm)^2/(2*\[Sigma]^2)))}, {sS, aS}, {z, L0, L1}, {t, -T, 0},
       MaxSteps -> Infinity]
    iniR = solS[[1]][[1]][[2]]
    Plot3D[Re[iniR[z, t]], {z, L0, L1}, {t, -T, 0}, PlotRange -> All]
    The following is the 2nd equation:
    Code (Text):

    solR = NDSolve[{D[sR[z, t], t] == (-(\[Gamma] - I*\[Beta]*z))*sR[z, t] - I*\[Omega]*aR[z, t],
        D[aR[z, t], z] + (1/c)*D[aR[z, t], t] == (-I)*\[Eta]*k*sR[z, t],
        sR[z, 0] == iniR[z, 0], aR[z, 0] == 0, aR[L0, t] == 0}, {sR, aR}, {z, L0, L1},
       {t, 0, T}, MaxSteps -> Infinity]
  13. Oct 28, 2011 #12
    Your code works. It's just there are numerical instabilities (or maybe the PDE itself is unstable) in the second PDE that get bad as t → T, and since you've set MaxSteps → Infinity, Mathematica keeps trying to refine the result. If you remove that option then the code runs fine.

    You should check the stability of the DE and if it is meant to be stable, then maybe playing with the Method option for NDSolve and the various stepsize and precision options will give you a reasonable solution of the the time range of interest...

    I've https://www.physicsforums.com/attachment.php?attachmentid=40440&stc=1&d=1319852925" to the post.

    Attached Files:

    Last edited by a moderator: Apr 26, 2017
  14. Oct 28, 2011 #13
    Simon: Thank you very much! It is good to know at least in principle it works in this way, though the solution is not good at all. And how do you know the instability happens at t-->T?
    Last edited by a moderator: Apr 26, 2017
  15. Oct 29, 2011 #14
    I'm just guessing from the plots of the numerical solution - run the code in the attached notebook. They get large and oscillatory along the z-axis as t gets large - not specifically as t → T. I'm a bit busy at the moment so I haven't actually looked closely at the PDE (it's also why it took me a while to get around to looking at your code).

    Your system is two first order, linear, coupled PDEs with complex coefficients. It should be amenable to analysis, but it's not really my area of expertise and off the top of my head, I'm not sure of the best way to go about it... Sorry.
  16. Oct 29, 2011 #15
    Simon: Yes, I looked at your attached notebook. Thank you very much! You already helped me a lot!
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook