Scalar PDE Coefficients in Function Form for the 'parabolic' pde solver

Hi, I want to use the pde parabolic solver for using a scalar parabolic pde in 2D. It is a'simple' Fokker-Planck equation and hence I have a first spatial derivative which the documentation suggests to put into the righthand side f. So I wrote f as function:

----------------------------------------- function f = fcoeff(p,t_mesh,u,t)

A=0.25; delta=0.25; omega=2*pi;

% Triangle point indices it1=t_mesh(1,:); it2=t_mesh(2,:); it3=t_mesh(3,:);

% Find centroids of triangles xpts=(p(1,it1)+p(1,it2)+p(1,it3))/3; ypts=(p(2,it1)+p(2,it2)+p(2,it3))/3;

x=xpts; y=ypts;

b1=-pi*A*sin(pi*(delta*sin(omega*t)*x.^2+(1-2*delta*sin(omega*t))*x)).*cos(pi*y); b2= pi*A*cos(pi*(delta*sin(omega*t)*x.^2+(1-2*delta*sin(omega*t))*x)).*sin(pi*y).*(2*delta*sin(omega*t)*x+1-2*delta*sin(omega*t)); size(u) [ux,uy] = pdegrad(p,t_mesh,u);

size(b1) size(ux)

f=b1.*ux+b2.*uy; --------------------------------------------------------------

and call the parabolic funkction via

u=parabolic(u0,tlist,'squareb1',p,e,t,c,@coefficient_a,@fcoeff,1);

It seems that the function fceoff does not get the matrices for u and the scalar for time t, as within the function they are 0 times 0 scalars. I guess I should use another notation?

[Bill Greene]

Hi,

I don't see anything obviously wrong with the code you posted.

Are you using the R2012b version or newer of MATLAB? Coefficients that are a function of u weren't supported before then.

If you want to post all of your code, I'll be happy to investigate further. If you do that, please use the {} Code format option so that I can extract the code in runnable form.

Bill

[mohammed aslam]

sir i need and coding of this eqation in matlab. it will be great help if u could do it..eqation 2.5 has to coded. i have tried every thing but i could'nt do it. i need an running program in mat. thank u..

Answer 1

[Andreas]

Hi Bill, I wrote this code on my Desktop PC with Matlab R2011a. As I compiled on a cluster with R2013a it worked. Thank you very much!!!!!!

I have a follow up question. I have to solve this scalar pde (second order parabolic with first derivatives ->Fokker Planck) for many initial densities. In the end, as this is a finite element code I think it's useful to compute the Mass and Stiffness matrix at every time step once and simply apply it to my different initial densities (and using my own implict Euler or something for time integration). Is this possible? I guess the elliptic solver assempde can put out these matrices, but if I got it right, with a coefficient depending on u or ux matlab wants to take the pdenonlin solver which cannot put out these matrices. Is there a possibility to get them anyway?

Answer 2

[Bill Greene]

Sorry, I'm not understanding your question.

If you want to compute K (or M or F) for a specific u-vector, you can call assema to do that for you. That is what the parabolic function is doing internally.

parabolic uses the ode15s function internally to integrate the equations in time. If you want to write your own fixed step size integrator, because you know a good step size or some other reason, you can probably use parabolic as a guide for how to do this. Note that parabolic does not reform the K or M matrices if they are not a function of time or the solution. Also, in my experience ode15s is fairly efficient so I don't think the time-step control is costing you that much.

Bill

Answer 3

[Andreas]

Sorry, if my question was a bit confused.

"Note that parabolic does not reform the K or M matrices if they are not a function of time or the solution." Yes, exactly! I have to solve my equation (what works now, thank you again) for many initial values. At the moment I do that simply via a for loop iterating over the initial densities.

    for i=1:N
        P(:,i)=parabolic(U0(:,i),....)
    end

Hence Matlab computes the same K and M matrices in every iteration. I would like to compute them just once and then use them to solve the equation for all the initial densities U0 at once. That is no problem in theory, but as my coeficients depend on u, I cannot use assempde but have to use pdenonlin (says so in the documentation, although my equation is linear), which does not put out the matrices K and M.

Is there a possibility to use assempde anyway?

Or can I even use parabolic directly for solving the equation for all the initial values at once? (via using a matrix or array containing the initial densities as initial value)?

Answer 4

[Bill Greene]

OK, so N is a very large number?

parabolic calculates a single time-dependent solution from a single initial solution vector.

You show above:

>u=parabolic(u0,tlist,'squareb1',p,e,t,c,@coefficient_a,@fcoeff,1);

You can use assempde to compute K, M, F, etc but internally it calls assema and assemb. You will need to call assema repeatedly to compute F since your f-coefficient is nonlinear but if you pass c=0 and a=0 it will skip the computation of K and M.

Doing something like this would probably allow you construct your own optimized, multi-solution form of parabolic.

Bill

Answer 5

[Andreas]

Thank you very much! One last question. I would need periodic boundary conditions. I have a two dimensional rectangle (e.g. [0,1]x[-8,8]) which is supposed to be periodic in x-direction. I am not sure if there is a way to implement this using parabolic. Is there a possibility to do that?

Answer 6

[Andreas]

ok, I guess there is no way then...

Answer 7

[mohammed aslam]

i am aslam. i just want to ask a running program of mat of the eqation that i have attached of eq 2.5. its a ballistic coefficient of reentry vehicle correction formula. i have tried all y best. partial derivative coding has to be done in mat. it will be great help if u could code it. also few thesis is mention in that page. thank u..