Mathematica: NDSolve with InterpolatingFunction as initial conditions

In summary, The questions arises because the user wants to use the solution of one PDE as an initial condition to solve another. They tried using NDSolve but it took a long time and caused issues with their computer's memory. Their questions are: 1. Can InterpolatingFunction be used as an initial condition in NDSolve? 2. Is there a way to bypass InterpolatingFunction and use a list of values as initial conditions? They have provided a simplified example of their problem.
  • #1
16
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
  • #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...
 
  • #3
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}
 
  • #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).
 
  • #5
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:
  • #6
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.
 
  • #7
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?
 
  • #8
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.
 
  • #9
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

  • NDSolve.nb
    15.2 KB · Views: 834
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!
 

What is Mathematica?

Mathematica is a computational software program used for mathematical and scientific calculations. It is commonly used by scientists, engineers, and mathematicians to solve complex equations and analyze data.

What is NDSolve in Mathematica?

NDSolve is a function in Mathematica that is used for solving ordinary differential equations (ODEs). It uses numerical methods to find approximate solutions to ODEs, which can then be used for further analysis.

What is an InterpolatingFunction in Mathematica?

An InterpolatingFunction is a type of data structure in Mathematica that represents a function by interpolating between a set of data points. It is often used as an initial condition in NDSolve to solve ODEs.

How do I use NDSolve with an InterpolatingFunction as initial conditions?

To use NDSolve with an InterpolatingFunction as initial conditions, you can first define the InterpolatingFunction and then use it as the initial condition in the NDSolve function. Alternatively, you can use the InterpolatingFunction directly in the NDSolve function by specifying it as the initial condition.

Can NDSolve with InterpolatingFunction as initial conditions handle boundary value problems?

Yes, NDSolve with InterpolatingFunction as initial conditions can handle boundary value problems (BVPs). BVPs involve finding a function that satisfies both the differential equation and certain boundary conditions. NDSolve uses numerical methods to find an approximate solution to the BVP.

Suggested for: Mathematica: NDSolve with InterpolatingFunction as initial conditions

Replies
2
Views
1K
Replies
4
Views
300
Replies
1
Views
70
Replies
1
Views
1K
Replies
1
Views
1K
Replies
2
Views
1K
Back
Top