pdepe extract intermediate values

Hi everyone,

I am solving a pdes using the pdepe function, but I have an issue.

Suppose you want to solve for u, and have your c, f, and s. However, s includes some function g(u)

Suppose that g is a sigmoid function that depends on the neighboring node.

For instance, g(i) = u(i) - sig(u(i-1) - threshold)

That way, the value of u at point i depends on neighboring nodes of previous time steps.

To overcome this issue, I would like to extract all my values of u at each time step into a function, calculate g and return back to the regular itteration. Basically calulate g at each time step, feed it into the giv equation, and proceed in solving.

Is there a way for me to extract all the solutions for u at each time step?

Thank you!

[Torsten]

Is there a way for me to extract all the solutions for u at each time step?

No, only the u-value(s) in the actual grid point x(i).

And the fact that you need the values of u in neighbouring nodes of previous time steps makes your problem a delay pde - even more complicated.

[Guy Elisha]

Hi, thanks for the response!

That is what I thought.

Also, sorry for being not clear, I dont need it in the previous time step but of neighboring nodes.

In my s term, I have s = (u * g) where g is a sigmoid function (g = tanh(u(i+1)-threshold)), so that if u ahead of the current node is greater than some threshold, g=1, and if u ahead of the current node is smaller than some threshold, g=-1.

Hence, s depends on u ahead of the current node (u(i+1))

To check that, I need basically all values of u at each time step.

Do you have a suggestion of how I can go around it?

Thanks again!

Answer 1

[Bill Greene]

I have written a PDE solver, pde1dm, that has some similarities to pdepe and accepts the same input syntax as pdepe. Most input that is acceptable to pdepe will also work in pde1dm.

However pde1dm has some additional capabilities. If you pass the option vectorized='on' to pde1dm, your pde function will be passed a vector of x-values spanning your complete spatial domain along with u and DuDx matrices at these x values. This allows your pde function to, for example, evaluate an integral involving the solution over the complete spatial domain. (Note that if you use this vectorized option, your pde function is also expected to return c, f, and s matrices corresponding to all values in the x vector.) You may be able to use this capability for your particular application.

You can find more information and download pde1dm here.

[Torsten]

@Bill Greene

Does your code change the (I guess) preassumed banded structure of the Jacobian in this case ?

[Bill Greene]

This vectorized option was added originally to improve performance, i.e. replace many calls to the user's pde function with a single call. So it's use (and implementation) has no effect on the structure of the equations.

[Torsten]

Usually, if the user has only access to the actual grid point, the Jacobian matrix for the integration has a banded structure and this is usually exploited in the solver setup.

But by having access to the complete solution vector, the user can include relations between the variables that destroy this banded structure, and the solver should be able to handle this case.

Or is the Jacobian in your code assumed to be full from the beginning ?

[Guy Elisha]

Interesting, thank you! I will check this out right now.

If I understand it correctly, this is exactly what I need. In my case, instead of solving an integral, I have a sigmoid that depends on u(i+1), but I think that it is a similar idea.

Thank you!

[Bill Greene]

This vectorized option assumes the jacobian is banded so it may only approximate the the "true" jacobian. There is a more advanced (and complicated) option for more general cases.

[Guy Elisha]

Hi @Bill Greene, I looked at your code, and I am still not sure how to extract the u values at any time step. Can you help me with that?

Thank you!

[Bill Greene]

The u values are available in the u-variable passed into your pde function but only at the single, specific time point requested by the ODE solver.

[Guy Elisha]

Thank you! Will try that.

I have two more questions:

1. when I run the heat conduction example, I get this error:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.011891e-17

Do you also encounter that?

2. Also, should I use pde1dM or pde1dm ?

[Bill Greene]

Post a complete mathematical description of the problem you are trying to solve along with your matlab code.

[Guy Elisha]

Before I even tried using your function on my code, I used one of your examples, for the heat conduction problem. This is where I got the error warning:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.011891e-17

Main code:

addpath('github_repo/')

%function heatConduction

% transient heat conduction in a rod

n=12; % number of nodes in x

L=1; % length of the bar in m

x = linspace(0,L,n);

t = linspace(0,2000,30); % number of time points for output of results

m=0; % rectangular Cartesian coordinate system

u = pde1dm(m, @heatpde,@heatic,@heatbc,x,t);

figure; plot(t, u(:,end)); grid on;

xlabel('Time'); ylabel('Temperature');

title('Temperature at right end as a function of time');

figure; plot(x, u(end,:)); grid on;

xlabel('x'); ylabel('Temperature');

title('Temperature along the length at final time');

%end

function [c,f,s] = heatpde(x,t,u,DuDx)

rho=8940; % material density

cp=390; % material specific heat

k=385; % material thermal conductivity

c = rho*cp;

f = k*DuDx;

s = 0;

end

function [pl,ql,pr,qr] = heatbc(xl,Tl,xr,Tr,t)

T0=100; % temperature at left end, degrees C

pl = Tl-T0;

ql = 0;

pr = 0;

qr = 1;

end

function u0 = heatic(x)

T0=100; % temperature at left end, degrees C

if x==0

u0=T0;

else

u0=0;

end

end

[Torsten]

Did you try to solve the system for g=0, e.g ?

[Guy Elisha]

@Torsten This was just an example. The system never actually yields something/zero, if this is what you are asking for.

I have a better explanation of the problem in my comment above (hand written notes). I am not giving so much information here just because it is a very convoluted problem that has a lot of layers into it. The main point is that I need all the values of u at every iteration so I can calculate g...

[Torsten]

I just asked if you tried to solve your problem without the sigmoidal function (or a function where you take its values in the current grid point).

I ask this because your problem looks so complicated that I would first try whether a simplified form is even solvable with pdepe/pde1dm.

[Guy Elisha]

You are right. So Previudly I had some pre-determined function for g (traveling sin wave), using pdepe and it works well. Now I am trying to add to this model by having g depend on the value of u instesd of a pre-assigned function.

Thank you by the way for your help!!

[Bill Greene]

Thank you for letting me know about this spurious warning message. It doesn't afftect the solution but it is clearly disconcerting. I uploaded a new version of the code that removes this message and, of course, you are free to download this. But, as I say, it doesn't afftect the results.

[Guy Elisha]

Hi @Bill Greene,

Unfortunetly, your function doesn't work to me... It is unclear what can be an input and what cannot.

For instance, why this is a valid input:

[u,vOde] = pde1dm(m,@pdeFunc,@icFunc,@bcFunc,x,t,@odeFunc,@odeIcFunc,xOde,opts);

and this is not:

[u,vOde] = pde1dm(m,@pdeFunc,@icFunc,@bcFunc,x,t,@odeFunc,@odeIcFunc,@FuncRand,xOde,opts);

where FuncRand is some function

Or, why is this not a valud input:

u = pde1dm(m, @heatpde,@heatic,@heatbc,x,t,@FuncRand);

Can you explain the inmut to pde1dm?

Thank you!

[Torsten]

Read Chapter 2 of the code documentation "pde1dM_manual.pdf": "Calling pde1dm"

It's a pdf-file in the folder "documents".

[Guy Elisha]

Thanks! Yes I saw that. Maybe my question was worded wrong.

In pdepe you can include variables and structs,

ex: sol = pdepe(0,@goveqns,@ICs,@BCs,xspan,tspan,a, b, c).

Is this an option in pde1dm?

[Torsten]

Pass extra arguments directly to the functions where they are needed, e.g.

sol = pdepe(0,@(x,t,u,dudx)goveqns(x,t,u,dudx,a,b,c),@ICs,@BCs,xspan,tspan)

That's the usual way in "modern" MATLAB and is also supported by pde1dm:

https://de.mathworks.com/help/optim/ug/passing-extra-parameters.html