MATLAB Solving Coupled Ingegro-Differential Equations with Maple/Matlab/Mathematica

[DrJekyll]

I've got a set of coupled Ingegro-Differential equations that I've been working on for a while. Here they are:

\frac{\partial b_x}{\partial t} = i x b_x + i u_x

\frac{\partial b_y}{\partial t} = i x b_y + i u_y

\frac{\partial b_z}{\partial t} = i x b_z + i u_z - \frac{3}{2 \omega} b_x

\frac{\partial u_x}{\partial t} = i x u_x + i b_x + \frac{2}{\omega} u_x -\frac{3}{2 \omega} \frac{\partial Q}{\partial x}

\frac{\partial u_y}{\partial t} = i x u_y + i b_y - \frac{1}{2 \omega} u_x -i Q

\frac{\partial u_z}{\partial t} = i x u_z + i b_z -\frac{i}{\kappa} Q

where Q is defined as follows:

Q = e^{-K x} \int_x^{\infty} \left( \frac{2 u_y(x')}{3} + \frac{2 i \omega u_x(x')}{9 K} \right) e^{K x'} dx' + e^{K x} \int_{-\infty}^x \left( \frac{2 u_y(x')}{3} - \frac{2 i \omega u_x(x')}{9 K} \right) e^{-K x'} dx'

Here I've included the dependence on x' in the integral to make it clear that the integral depends on ux and uy.

and K = \frac{2}{3} \omega \sqrt{1+\frac{1}{\kappa^2}}

\kappa and \omega are parameters and if need be can be set to 0.0001 and 0.01 respectively.

To date I've written a code in Fortran 90 to solve these equations using a 4rth order Runge-Kutta method along with an adaptive Simpson integration. The code runs and spits out results but I question the validity. Thus, I would like to try and get Maple/Matlab/Mathematica to run through at least one time step to compare too (mainly the result of the integration).

I'm confident that Maple just can't do it, and I'm unsure of the exact syntax to make Mathematica work. As for Matlab, I've been messing around with the pdepe function with an integration call plopped into it but I've had little sucess.

If anyone can help with this it would be much appreciated.

 

[kosovtsov]

If I understand your notation correctly

u=u(t,x,y,z), b=b(t,x,y,z)

and, e.g.,

u_x=\frac{\partial u}{\partial x}

then your system has only trivial solution

b(t,x,y,z) = _F1(t),u(t,x,y,z) = _F2(t)

 

[DrJekyll]

My notation may be a bit confusing. b_x represents the x component of b not the partial derivative. Also b=b(x,t), the y and z directions can be safely ignored.

 

[kosovtsov]

I believe that it would be better to convert your Integro-Differential system into purely (linear) PDEs system replacing the last condition on Q by the following PDE

-9K\frac{\partial^2 Q}{\partial x^2}+9K^3Q+12K^2u_y-4i\,\omega\frac{\partial u_x}{\partial x}=0

with the boundary conditions Q|_{x=\pm \infty}=0.

I think that there is a chance to solve your system exactly (symbolicaly).

 

[DrJekyll]

I've just recently tried this as well. One difference though,

\frac{\partial^2 Q}{\partial x^2} - K^2 Q = \frac{4 \partial u_y}{3 \partial x} + \frac{4 i \omega u_x}{9}

(without actually checking, they look mostly the same, the above is definitely the correct one, any differences would most likely be my fault in the original post)

It should be noted that the above is a constraint equation.

I still can't get Maple or Mathematica to solve with the constraint equation in. Matlab does produce a result which I'm working on now, resolution might be a problem.

 

[DrJekyll]

Ok, so resurrecting this thread since I'm back to trying to solve using Matab.

Here are the equations I'm trying to solve again, with the last one being a constraint equation. So far I've been using a 4rth order Runge-Kutta scheme to solve the first 6 equations with some initial conditions for 1 time step. Then I've used another RK4 scheme to solve the last equation (after breaking it into 2 first order ODE's) using the values from the previous RK4 step.

\frac{\partial b_x}{\partial t} = i x b_x + i u_x

\frac{\partial b_y}{\partial t} = i x b_y + i u_y

\frac{\partial b_z}{\partial t} = i x b_z + i u_z - \frac{3}{2 \omega} b_x

\frac{\partial u_x}{\partial t} = i x u_x + i b_x + \frac{2}{\omega} u_x -\frac{3}{2 \omega} \frac{\partial Q}{\partial x}

\frac{\partial u_y}{\partial t} = i x u_y + i b_y - \frac{1}{2 \omega} u_x -i Q

\frac{\partial u_z}{\partial t} = i x u_z + i b_z -\frac{i}{\kappa} Q

<br /> \frac{\partial^2 Q}{\partial x^2} - K^2 Q = \frac{4 \partial u_y}{3 \partial x} + \frac{4 i \omega u_x}{9} <br />

with
<br /> K = \frac{2}{3} \omega \sqrt{1+\frac{1}{\kappa^2}} <br />

So, as I said I'm trying to use Matlab to solve this now. Previously I tried using Matlabs pdepe function to solve it and I got a result which I don't believe. Now I would like to use the ODE functions to try and solve it. My guess is that I will need to use something similar to my code to solve it, in that I'll need to solve the first 6 equations for 1 time step, then solve the constraint equation.

The problem is that I have no idea how to do this in Matlab, I have next to no experience with it so any help would be appreciated.