Mathematica: NDSolve with InterpolatingFunction as initial conditions

[crazybird]

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?

 

[Simon_Tyler]

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...

 

[Bill Simpson]

2. Is there any way to ...make a list of values out of InterpolatingFunction...

Here is an example interpolating function and a table created from it.

In[1]:=f=Interpolation[Evaluate[Table[{x,1/x},{x,1,5}]]];
Table[f,{i,2,6}]

Out[1]= {1/2, 1/3, 1/4, 1/5, 2/15}

 

[FunkyDwarf]

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).

 

[crazybird]

Hi, All, here is my example:
c = 29.979;
T = 500;
tm = -250;
σ = 100;
β=0.133427;
γ=10^-5;
w=1.32949*10^-8;
h=1539.06;
k=40895;

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σ²), aS[L0,t]==E^((t-tm)²/2σ²),{sS,aS},{z,L0,L1},{t,-T,0}];

iniR=solS[[1]][[1]][[2]];

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}]

 

[crazybird]

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).

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.

 

[crazybird]

Here is an example interpolating function and a table created from it.

In[1]:=f=Interpolation[Evaluate[Table[{x,1/x},{x,1,5}]]];
Table[f,{i,2,6}]

Out[1]= {1/2, 1/3, 1/4, 1/5, 2/15}

Bill: Yes, one can make a list like this. But the question is how to make it into an initial condition of NDSolve?

 

[crazybird]

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...

Simon: Example is given in 5th floor.

 

[jackmell]

Hi, All, here is my example:
c = 29.979;
T = 500;
tm = -250;
σ = 100;
β=0.133427;
γ=10^-5;
w=1.32949*10^-8;
h=1539.06;
k=40895;

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σ²), aS[L0,t]==E^((t-tm)²/2σ²),{sS,aS},{z,L0,L1},{t,-T,0}];

iniR=solS[[1]][[1]][[2]];

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}]

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.

 

[crazybird]

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.

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.

 

[crazybird]

Working Mathematica Note:

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

\[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:

(**********************)
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]

 

[Simon_Tyler]

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.

 

[crazybird]

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.

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?

 

[Simon_Tyler]

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.

 

[crazybird]

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.

Simon: Yes, I looked at your attached notebook. Thank you very much! You already helped me a lot!