How to control times at which PDE Toolbox solves a parabolic equation?

I am solving a 1D groundwater flow equation with time-dependent forcing using PDE Toolbox. I managed to input the time-dependent forcing (which is a vector of data values at specific times) into the PDE Toolbox using the following function:

   function f = fcoeff(t)

   W = importdata('Rainfall_daily.mat'); % the forcing 
   nt=length(W);
   tW = 0:nt-1; % values of time at which forcing data is available
   W_intrp = interp1(tW,W,t,'linear'); % interpolated values of forcing at 
   times used by toolbox
   disp([t, W_intrp])
   f = W_intrp;

However, the PDE toolbox misses some data points in the vector 'W' because it sometimes takes too large time steps. This does not yield correct results at the time values I ask it to plot the results at. That is why I want to force the PDE toolbox to solve the parabolic PDE at time values that I specify. This will also help me avoid using interp1 in the above function.

Thanks, Abrar

Answer 1

[Bill Greene]

I'm afraid I wasn't able to run your example. I believe that your fcoeff function has an error in it.

I did take a look at the Rainfall_daily.mat file. If I am right in assuming that the first column is some kind of time measure and the second is the function value, it looks like the function is extremely discontinuous.

The ODE solver used by the PDE Toolbox parabolic function has an option, MaxStep, that allows you to specify the maximum step size used during the solution. It might be possible to improve the solution you are getting by setting this to a fraction of the delta_time of your data.

Unfortunately, there is no simple way to do this without editing the PDE Toolbox parabolic function. If you edit this function (I suggest making a copy of the original first), you will see where the ode15s function is being called and where the other options to ode15s are being set. Adding the MaxStep option should be straightforward. I think this change will cause the parabolic function to execute much slower so I suggest you remove this change when you are through with this particular example.

Good luck!

Bill

Answer 2

[Bill Greene]

>This does not yield correct results at the time values I ask it to plot the results

So, why don't you request results at the time of each point in your W-vector?

You could also experiment with setting smaller values of the rtol and atol optional input arguments to the parabolic function.

Regardless, you will still need the interp1() function to perform interpolation because the function ode15s, the ODE solver used by parabolic(), will call your function at arbitrary times between the initial and final times.

Bill

[Abrar Habib]

I do request time at the specific times that I have data for and it does output results at the times I asked for, but because I have a disp() function in my fcoeff function I can see that it skips some values in my W vector and this explains the results I get out of the PDE toolbox where I should see, say a peak in the plot which is not displayed in the results. I assume it somehow interpolates the results to produce them at where I want them exactly after performing computation.

I haven't tried playing with rtol and atol but I will give it a shot.

If I manage to find a way to force the PDE toolbox to perform computations at the times for which I have forcing data (i.e. my W vector) then I can simply take the values from the W vector without interpolating which I would prefer because the interpolation assumes a certain distribution of the forcing over a time step(in my case I assume it is linear) which is not true and hence I introduce an assumption that I hope I can avoid.

Thank you for the reply, Abrar

[Bill Greene]

So if I understand you right, you are including every time in your W vector in the tlist argument to parabolic, but parabolic is not calling your f-coefficient evaluator at those times?

I would need to investigate this further. If you can post your complete example, (using the paper-clip icon in the editor dialog below), I'll take a look.

Bill

[Abrar Habib]

Thanks for your reply Bill and sorry for the late response.

I have been playing around with the atol and rtol values and this helps by reducing the time step values (obviously)and hence taking into account most rainfall events, however as I mentioned earlier this introduces inaccuracies to my system where, due to interpolation in the fcoeff function, additional rainfall is accounted for. So after days of playing around with the fcoeff function to try and reduce this inaccuracy, I think I am back to square one where I think the solution to this problem is to dictate the times at which PDE Tool performs computation (if this is possible).

I have attached both the PDE model and a function, in addition to the fcoeff pasted in the first post, that is used by the model, and the rainfall data file I use.

Thanks for your effort, Abrar