Mathematica Mathematica: NDSolve with InterpolatingFunction as initial conditions

AI Thread Summary
Using InterpolatingFunction as an initial condition in NDSolve is feasible, but it requires pre-evaluating the function to ensure that the inputs are numerical values. Users have suggested creating a list of values from the InterpolatingFunction to use as initial conditions, which can be implemented using loops or tables. However, it is crucial to ensure that the initial conditions and boundary conditions are correctly specified for both functions involved in the PDEs. Numerical instabilities may arise in the second PDE, particularly as time approaches the upper limit, which can be mitigated by adjusting the method and precision options in NDSolve. Overall, careful handling of initial conditions and stability analysis is essential for successful implementation.
crazybird
Messages
16
Reaction score
0
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?
 
Physics news on Phys.org
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...
 
crazybird said:
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}
 
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).
 
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/2σ2), aS[L0,t]==E^((t-tm)2/2σ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}]
 
Last edited:
FunkyDwarf said:
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.
 
Bill Simpson said:
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?
 
Simon_Tyler said:
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.
 
crazybird said:
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/2σ2), aS[L0,t]==E^((t-tm)2/2σ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.
 
  • #10
jackmell said:
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.
 
  • #11
Working Mathematica Note:

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

Code:
\[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:
(**********************)
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]
 
  • #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.
 

Attachments

Last edited by a moderator:
  • #13
Simon_Tyler said:
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?
 
Last edited by a moderator:
  • #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.
 
  • #15
Simon_Tyler said:
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!
 

Similar threads

Replies
1
Views
2K
Replies
4
Views
3K
Replies
6
Views
2K
Replies
8
Views
5K
Replies
3
Views
5K
Replies
5
Views
3K
Back
Top