How to solve below PDE based model using finite difference scheme.

I want to solve a PDE system which is given by

with initial conditions

I want to apply finite difference method to solve the above system for f and g.

We discretize f(t,a,x) as f(i,j,k) and same goes for function "g" as well. If we apply forward finite difference scheme, then both the partial derivative becomes

I know if I start with intial condition then we can solve this. From initial condition we have f(1,j,k)=f(0,a,x) and g(1,j,k)=g(0,a,x) in our hand, putting both in main equation we can get f(2,j,k) and g(2,j,k), but I don't how to handle those integral terms given in both equations . Please Can anyone with PDE background help me.

Any idea about this will help me a lot.

[Torsten]

And what is the a-x-domain over which you want to solve your equations ?

[Nikhil]

@Torsten and . If possible can you give some hints.

[Bill Greene]

Then why are the integrals from 0 to 40 in your equation definition?

[Nikhil]

@Bill Greene Sorry It was a typing mistake. I have updated the x-domain now to [0,40]

[Torsten]

May I ask where these equations stem from ?

Is it an artificial problem or does it have a physical background ?

[Bill Greene]

One rather fundamental problem you have is that your initial conditions are inconsistent.

For example, just consider the point t=0, a=0, x=0.

[Nikhil]

@Torsten these type of equations are mostly used in population dynamics.

[Nikhil]

@Bill Greene I agree with what you are saying but the thing is we are going to solve these equations in MATLAB. so at t=x=a=0 we get g(0,0,0)=0.0009 and according to me MATLAB will take it as ZERO itself . I could have taken to answer your query. But then a*e^{-0.08a} wil not decrease as rapidly as the function which I had taken.

[Nikhil]

@Torsten and @Bill Greene do you know how to deal with the integral part involving in these equations.

Answer 1

[Torsten]

Take a grid (ai,xj) over your rectangular region with

0 = a1 < a2 < ... < an = 60

0 = x1 < x2 < ... < xm = 40

Approximate the integrals using the trapezoidal rule by using the values in the grid points.

Then you get a system of 2*(n*m) ordinary differential equations for f and g in the grid points that can be solved using ode15s.

I couldn't classify your problem into a certain category (it's not a PDE as you claim). So I don't know if the method I suggested will work. If you have literature about numerical methods to solve such systems (this is the reason I asked for the application), you should study it first before wildly beginning to produce MATLAB code.

[Bill Greene]

Yes, I think that is a very good approach. The trapz function should make it straightforward to calculate the four required integrals.

[Nikhil]

@Torsten I also tried using trapezoidal rule but as you can see in the second integral , here the integrand is function of both x and y and both variable runs from 0 to 40 and we are integrating the function with respect y only. How to use trapezoidal in such cases.

[Torsten]

I don't understand the problem.

You are given f(t,ai,xj) for 1<= i<= n and 1<=j<=m from the integrator. Now you want to evaluate the integral for a certain index pair (I,J) at (t,aI,xJ) with 1<=I<=n and 1<=J<=m

Thus it's approximately

exp(-(xJ-x(1))^2)*f(t,aI,x(1))/sqrt(pi) * dx/2 + 
(exp(-(xJ-x(2))^2)*f(t,aI,x(2))/sqrt(pi) + exp(-(xJ-x(3))^2)*f(t,aI,x(3))/sqrt(pi) + ... + exp(-(xJ-x(m-1))^2)*f(t,aI,x(m-1))/sqrt(pi)) * dx +
exp(-(xJ-x(m))^2)*f(t,aI,x(m))/sqrt(pi) * dx/2