Help with Maple/Matlab/Mathematica

[Biftheunderstudy]

I have a set of 2 coupled differential vector equations. When I expand them all out they form a set of 6 coupled 1st order DE's (they are complex).

This is a bit of a rehash of an earlier post from way back when.
Long story short, I've been working on solving these equations numerically using Fortran90 and I've been getting funny results. To test my solution I've been toying with Maple/Matlab/Mathematica to see if it can solve these equations, after posting a few times on their respective websites I've gotten nowhere. Only the Maple forums actually responded and basically told me that Maple was ill-suited to solve this problem.

Here is exactly what I've put into Mathematica and it yells at me for encountering non-numerical value for a derivative.

If Matlab can solve this quicker/easier I'm open to suggestions. Help is much appreciated.

Now I'll explain what everything is:

ux...uz are velocity components, bx..bz are magnetic field components

Q is a pressure, calculated from a divergence equation. It comes out as a solution to a Helmholtz equation and takes the form of an integral.

k is the constant for the green's function Q'' - k^2 Q = F

kappa and omega are parameters to be varied

x (really should be called r or something to be less confusing) is the spatial coordinate

sigma is just there for the boundary conditions

The boundary conditions are 0 everywhere, with a Gaussian in bx and a corresponding wiggle in by to make the divergence 0.

k := 2/3 * \[Omega] * Sqrt[1 + 1/\[Kappa]^2]

IC[x_] := E^(-(x^2)/2/\[Sigma]^2)

\[Sigma] := 10./(2.*Sqrt[2.*Log[2]])

dQ[x_] := Integrate[(2*k*uy[t, y]/3 + 2*I*\[Omega]*ux[t, y]/9)*
E^(-k*y), y]*E^(k*x) +
Integrate[(-2*k*uy[t, y]/3 + 2*I*\[Omega]*ux[t, y]/9)*E^(k*y), y]*
E^(-k*x) + 4*uy[t, x]/3

Q[x_] := Integrate[(2*uy[t, y]/3 + 2*I*\[Omega]*ux[t, y]/9/k)*
E^(-k*y), {y, x, 50}]*E^(k*x) +
Integrate[(2*uy[t, y]/3 - 2*I*\[Omega]*ux[t, y]/9/k)*
E^(k*y), {y, -50, x}]*E^(-k*x)

\[Kappa] := 0.01
\[Omega] := 0.01

sol = NDSolve[
{D[bx[t, x], t] == I*x*bx[t, x] + I*ux[t, x],
D[by[t, x], t] == I*x*by[t, x] + I*uy[t, x] - 3*bx[t, x]/\[Omega]/2,
D[bz[t, x], t] == I*x*bz[t, x] + I*uz[t, x],
D[ux[t, x], t] ==
I*x*ux[t, x] + I*bx[t, x] + 2*uy[t, x]/\[Omega] - 3*dQ/2/\[Omega],
D[uy[t, x], t] ==
I*x*uy[t, x] + I*by[t, x] - ux[t, x]/2/\[Omega] - I*Q,
D[uz[t, x], t] == I*x*uz[t, x] + I*bz[t, x] - I*Q/\[Kappa],
bx[0, x] == E^(-(x^2)/2/\[Sigma]^2),
by[0, x] == 3*I*D[E^(-(x^2)/2/\[Sigma]^2), x]/2/\[Omega],
bz[0, x] == ux[0, x] == uy[0, x] == uz[0, x] == 0},
{bx, by, bz, ux, uy, uz}, {t, 0, 10}, {x, -50, 50}]

 

[CycloWizard]

When you say they're complex, you mean that they include imaginary components or that they're really complicated? I assume it's the former. I'm not terribly familiar with Mathematica for solving numerical problems, but I know that this problem can be solved in MATLAB using the 'pdepe' function. Let me know if you need help figuring out how to set it up in more detail - it's always a pain to set up MATLAB's DE solvers.

 

[Biftheunderstudy]

Its both complicated and complex

I am completely new to Matlab, so any and all help will be appreciated on that front.

 

[schenley101]

I believe that it is possible to do with matlab using the symbolic toolbox. It can do ODEs that don't have numerical answers. The problem with the symbolic tool box is that it is not always 100% correct. Do you already have the analytical solution to the problem?

 

[schenley101]

I would also check out matlab central to see if anyone else has solved a similar problem

 

[CP5670]

I think there is a problem with how you are using Q and dQ. You have defined them to be functions of x, but you use them in NDSolve like constants.

Mathematica and Matlab can both do problems like this. As a general rule, I use Mathematica for doing quick-and-dirty calculations, and Matlab for production-level programs. Mathematica's symbolic capabilities allow you to be sloppy with a lot of things that Matlab complains about and spend far less time programming, but that can also make your computations slow. I never used Maple much and can't comment on that.

To test my solution I've been toying with Maple/Matlab/Mathematica to see if it can solve these equations, after posting a few times on their respective websites I've gotten nowhere. Only the Maple forums actually responded and basically told me that Maple was ill-suited to solve this problem.

I once posted on the Mathematica forum about something, and I actually got a response after 7 years.

 

[CycloWizard]

Really quick since I'm busily working today (for once)... This solves a diffusion equation for two species with no source term. I'll try to take a more detailed look later, but this should get you pointed in the right direction anyway.

%Main function call to pdepe:
w=pdepe(0,@pdefun,@ICFun,@BCFun,z,t);

function [c,f,s] = pdefun(x,t,u,DuDx)
%Compute the values of the PDE components c, f, and s, where
% c is the coefficient of the time derivative
% f is the coefficient of the Laplacian
% s is the source term
c=[1;1];
f=[DuDx(1);0];
s=[0;0];

function w=ICFun(x)
%gives initial conditions
w=[1;1];

function [pl,ql,pr,qr] = BCFun(xl,ul,xr,ur,t)
%Return the boundary conditions for input to the pde solver pdepe.
pl=[1;0];%constant concentration of species 1 on left side
ql=[0;1];%zero flux of species 2 at left side
pr=[0;1];%constant concentration of species 2 on right side
qr=[1;0];%zero flux of species 1 at right side

 

[canis]

..

 

[johnnydrinkawater]

I've always used the symbolic math toolbox to solve diff eq's in the past as a stepping stone and separate m-file. I could never figure out how to get the numerical diff functions to work right for me.

Copy+paste symbolic solutions into the numerical m-file.

 

[CycloWizard]

A problem like this likely does not have an analytical solution guys. Most problems don't, at least not the interesting ones. If they do, they've already been solved.

 

[Biftheunderstudy]

Cyclowizard is right, no analytical solution.
I'll get my hands on a copy of Matlab today and start giving it a whirl. Never used it before so it should be a journey....

Just to make sure you didn't miss a key point though. The variable Q is defined as an integral, and it has an exponential which means it can not be differentiated away. This makes this problem very difficult. That integral has been the source of so much frustration as it takes the longest time to compute in my code.

 

[CycloWizard]

Cyclowizard is right, no analytical solution.
Just to make sure you didn't miss a key point though. The variable Q is defined as an integral, and it has an exponential which means it can not be differentiated away. This makes this problem very difficult. That integral has been the source of so much frustration as it takes the longest time to compute in my code.

That's not a major obstacle in MATLAB: just solve for Q numerically using the quad function (or trapz, if it's more convenient). This can be integrated in a straightforward way in your script. Like I said, let me know if you get stuck and I'll take a closer look.

 

[Biftheunderstudy]

Alright, I've been messing around with Matlab for a bit and I've put most of my problem into the form of a system of diffusion equations. Here is what I've put into the M-file.

function LinearEquationSolving
x=linspace(-50,50,100);
t=linspace(0,5,500);

%Main function call to pdepe:
sol=pdepe(0,@pdefun,@ICFun,@BCFun,x,t);
B1 = sol,:,1);
B2 = sol,:,2);
B3 = sol,:,3);
U1 = sol,:,4);
U2 = sol,:,5);
U3 = sol,:,6);

function [c,f,s] = pdefun(x,t,u,DuDx)
%Compute the values of the PDE components c, f, and s, where
% c is the coefficient of the time derivative
% f is the coefficient of the Laplacian
% s is the source term
c=[1;1;1;1;1;1];
f=[0;0;0;0;0;0];

omega = 0.01;
kappa = 0.00001;
K = 2*omega*sqrt(1+1/kappa^2)/3;

F1 = 1i*x*u(1) + 1i*u(4);
F2 = 1i*x*u(2) + 1i*u(5) - 1.5*u(4)/omega;
F3 = 1i*x*u(3) + 1i*u(6);

dQ = exp(K*x)*quad(@dQintl,x,100) + exp(-K*x)*quad(@dQintr,0,x)+4*u(5)/3;
Q = exp(K*x)*quad(Qintl,x,100) + exp(-K*x)*quad(Qintr,0,x);

F4 = 1i*x*u(4) + 1i*u(1) + 2*u(5)/omega - 1.5*dQ;
F5 = 1i*x*u(5) + 1i*u(2) -0.5*u(4)/omega - 1i*Q;
F6 = 1i*x*u(6) + 1i*u(3) - 1i*Q/kappa;

s=[F1;F2;F3;F4;F5;F6];

EDIT:

function answer=dQintl(x,u)
omega = 0.01;
kappa = 0.00001;
K = 2*omega*sqrt(1+1/kappa^2)/3;
answer = (2i*omega*u(1)/9 - 2*u(5)*K/3)*exp(K*x);
function answer=dQintr(x,u)
omega = 0.01;
kappa = 0.00001;
K = 2*omega*sqrt(1+1/kappa^2)/3;
answer = (2i*omega*u(1)/9 + 2*u(5)*K/3)*exp(-K*x);
function answer = Qintl(x,u)
omega = 0.01;
kappa = 0.00001;
K = 2*omega*sqrt(1+1/kappa^2)/3;
answer = (2*u(5)/3 - 2i*omega*u(4)/(9*K))*exp(K*x);
function answer = Qintr(x,u)
omega = 0.01;
kappa = 0.00001;
K = 2*omega*sqrt(1+1/kappa^2)/3;
answer = (2*u(5)/3 + 2i*omega*u(4)/(9*K))*exp(-K*x);

\End EDIT

function w=ICFun(x)
%gives initial conditions
omega = 0.01;
kappa = 0.00001;
sig = 10./(2.*sqrt(2.*log(2.)));
bx0 = exp(-(x^2)/(2*sig^2));
by0 = 3i*x*exp(-(x^2)/(2*sig^2))/(2*sig*omega);
w=[bx0;by0;0;0;0;0];

function [pl,ql,pr,qr] = BCFun(xl,ul,xr,ur,t)
%Return the boundary conditions for input to the pde solver pdepe.
pl=[0;0;0;0;0;0];%constant concentration of species on left side
ql=[0;0;0;0;0;0];%zero flux of species 2 at left side
pr=[0;0;0;0;0;0];%constant concentration of species 2 on right side
qr=[0;0;0;0;0;0];%zero flux of species 1 at right side


I'm still unsure of how to make the boundary conditions correct, I want them to be just 0 everywhere on the boundary.

As for the rest of it, there seems to be a problem with how I'm passing the quad function. I think this arises because the integral needs to know ux(x') and uy(x') and I'm not sure how to make Matlab understand this. (x' is the integration variable)

Should I be using xmesh) or something like that when I want to use the coordinate variable x?

Here is the list of error messages I get when I evaluate it:

??? Undefined function or variable 'x'.

Error in ==> LinearEquationSolving>dQint at 43
if y<x

Error in ==> quad at 77
y = f(x, varargin{:});

Error in ==> LinearEquationSolving>pdefun at 30
dQ = exp(K*x)*quad(@dQint,-50,x) + exp(-K*x)*quad(@dQint,x,50)+4*u(5)/3;

Error in ==> pdepe at 259
c = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in ==> LinearEquationSolving at 6
sol=pdepe(0,@pdefun,@ICFun,@BCFun,x,t);

 

[drinkmorejava]

You could probably try ODE45 (Matlab). I'm a bit busy right now, so I can't help too much, but that will solve as many coupled ODEs as you need. If it's the first time you've used it though, it'll probably be a bitch to figure out.


LOL, nm. I probably could have read the thread and seen that you have PDEs.

 

[Biftheunderstudy]

I have 2 main problems left unsolved.

1) I don't have a firm grasp on how the boundary conditions should be input. Basically, all of the variables should smoothly go to zero at the left boundary (x=0) and the right boundary (x=100). I've taken this to mean that pl,ql,pr,qr are all =0 but I'm pretty sure this is wrong.

2) The integral is really bugging me. The function quad has 3 main arguments, it wants the function handle that you want to integrate and the upper and lower bounds.

a)The first argument is the function handle, which needs to take a parameter x and return a value. This is where I run into some trouble. If the integration variable is x' the integral depends on ux(x') and uy(x'). This makes it difficult since I can't define a function that would only depend on x, I would need to make it depend on ux and uy as well. Something like:

function answer = Qint(x,u)

I think the solution might lie somewhere with nested functions but I know absolutely nothing about these.

b) The bounds of the integrals are from 0 to x and from x to 100. This is a truncation of -inf to inf which then went to -inf to x and x to inf (I've truncated infinity to 0 to 100)

 

[CycloWizard]

I have 2 main problems left unsolved.

1) I don't have a firm grasp on how the boundary conditions should be input. Basically, all of the variables should smoothly go to zero at the left boundary (x=0) and the right boundary (x=100). I've taken this to mean that pl,ql,pr,qr are all =0 but I'm pretty sure this is wrong.

I would recommend looking at the help file for the pdepe function - it will give you the exact form of the BCs. IIRC, the L refers to left boundary (x=0) and R is the right (x=100); p represents the Dirichlet component of the BC; and q represents the coefficient function of the spatial gradient (their product combined would be the flux).

2) The integral is really bugging me. The function quad has 3 main arguments, it wants the function handle that you want to integrate and the upper and lower bounds.

a)The first argument is the function handle, which needs to take a parameter x and return a value. This is where I run into some trouble. If the integration variable is x' the integral depends on ux(x') and uy(x'). This makes it difficult since I can't define a function that would only depend on x, I would need to make it depend on ux and uy as well. Something like:

function answer = Qint(x,u)

I think the solution might lie somewhere with nested functions but I know absolutely nothing about these.

b) The bounds of the integrals are from 0 to x and from x to 100. This is a truncation of -inf to inf which then went to -inf to x and x to inf (I've truncated infinity to 0 to 100)

I think you'll have to use nested functions. Generally, these have the form
Integral_Value = quad(@(x) Qint(x,u),lim1,lim2)
thus, @(x) is the variable on which the quad would operate. I haven't actually had any reason to use quad before, so I'm not sure if this is the exact syntax. Post back if you can't get it working and I'll take a look at it later tonight when I'm not quite as busy.

 

[Biftheunderstudy]

Using the nested function seemed to make it do something, though it just yells that everything is NaN.

As for the BC's, I have Dirichlet BC's so q=0, and, if my digging through tons of webpages on this, then p = u1(1..6). My BCfun should look like this then:

function [pl,ql,pr,qr] = BCFun(xl,ul,xr,ur,t)
&#37;Return the boundary conditions for input to the pde solver pdepe.
pl=[ul(1);ul(2);ul(3);ul(4);ul(5);ul(6)];%constant concentration of species on left side
ql=[0;0;0;0;0;0];%zero flux of species 2 at left side
pr=[ur(1);ur(2);ur(3);ur(4);ur(5);ur(6)];%constant concentration of species 2 on right side
qr=[0;0;0;0;0;0];%zero flux of species 1 at right side


The NaN's are most likely caused by the integral acting funny and I'm pretty sure I know where the issue is.
Basically I need to do this integral:

exp(K*x)*int_0^x (uy(x',t) + i*ux(x',t))*exp(-K*x') + exp(-K*x)*int_x^100 (uy(x',t) -i*ux(x',t))*exp(K*x')

having x in the bounds of integration might be causing problems. Also, I'm not sure if simply passing the function u will make it integrate it properly.

 

[CycloWizard]

Oops, missed your reply here. I think your boundary conditions are fine. I'd recommend looking here for the documentation on 'quad':
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/quad.html

"q = quad(fun,a,b) tries to approximate the integral of function fun from a to b to within an error of 1e-6 using recursive adaptive Simpson quadrature. fun is a function handle. See Function Handles in the MATLAB Programming documentation for more information. Limits a and b must be finite. The function y = fun(x) should accept a vector argument x and return a vector result y, the integrand evaluated at each element of x."

Having x in the function itself shouldn't be a problem, as long as the limits a and b are finite constants. I think you said you were using 0 and 100 for a and b, so that shouldn't be an issue. You might try setting up a simpler trial function to ensure you have the syntax right for the parametrization of the quad function, as that can be a bit tricky. You might also need to use the quadv function (vectorized quadrature), but it's unclear to me exactly what you're passing to this function at this point.

 

[Biftheunderstudy]

The original integral goes like exp(-K|x-x'|) from -infinity to infinity.
After splitting the absolute value and truncating it becomes:
exp(-K (x-x') from 0 to x and exp( K (x-x')) from x to 100

Does the quad function know if I pass it x as a bound that it should evaluate that at a particular scalar x? Same with the exp(-K*x) out front? My thoughts are that it will interpret them as vectors and compute the entire Q vector all at once, but I consider this unlikely.

To fix this, I would want to make a for loop like:

for n=0:end
xd = x(n)
Q(n) = exp(-K*xd)*quad(@(x)Qintl(x,u),xd,100) + exp(K*xd)*quad(@(x)Qintr(x,u),0,xd)
end
(I think most of this loop is correct except the n=0:end part. How do I tell it to loop till the end of the array?)

Also, I suspect that when I pass the quad function u, it doesn't know that it should be evaluating u(x') and I'm unsure as to how to do this.

P.S.
Adaptive simpson quadrature is actually the same algorithm that my code has been using, except parallel.

 

[canis]

..

 

[CycloWizard]

The original integral goes like exp(-K|x-x'|) from -infinity to infinity.
After splitting the absolute value and truncating it becomes:
exp(-K (x-x') from 0 to x and exp( K (x-x')) from x to 100

Does the quad function know if I pass it x as a bound that it should evaluate that at a particular scalar x? Same with the exp(-K*x) out front? My thoughts are that it will interpret them as vectors and compute the entire Q vector all at once, but I consider this unlikely.

To fix this, I would want to make a for loop like:

for n=0:end
xd = x(n)
Q(n) = exp(-K*xd)*quad(@(x)Qintl(x,u),xd,100) + exp(K*xd)*quad(@(x)Qintr(x,u),0,xd)
end
(I think most of this loop is correct except the n=0:end part. How do I tell it to loop till the end of the array?)

Also, I suspect that when I pass the quad function u, it doesn't know that it should be evaluating u(x') and I'm unsure as to how to do this.

P.S.
Adaptive simpson quadrature is actually the same algorithm that my code has been using, except parallel.

The loop would be
for n=0:length(array)
Remember that MATLAB uses indexing starting at 1 (rather than at 0 like C does). You can also simply use
for n=indices
where "indices" is an array containing all of your indices.

As far as whether something is evaluated as a scalar or vector, it depends on the notation you use. If you need to use element-wise multiplication/division/exponents, you need to use .*, ./, or .^ rather than *, /, or ^. It should know what to do with the input you gave above as long as xd is a scalar and u is a vector the same length as x (what is u again? ).

 

[Biftheunderstudy]

Thanks, the ./, .* etc was one of the key parts. With a bit more work I realized that I can use a constraint equation instead of the integral which is of the form:

Q'' - K^2 Q -F = 0

where F is a function of ux and uy.

With this Maple and Mathematica still fail, Matlab however was able to solve this with pdepe. I'm still analyzing the solution, but it looks wrong

Any advice on how to write a simple finite difference code? Mainly, do I have to do anything differently to deal with a system of equations?
example:

bx(x,t+dt) = dt*(i*x*bx(x,t) + i*ux(x,t)) + bx(x,t)

 

[CycloWizard]

It's very easy to write finite difference code in Matlab, but writing an efficient one is a bit trickier. The key is to ensure everything is vectorized as much as possible. So, rather than doing element-wise computation for each spatial node at each time step, all of the spatial nodes should be computed in the same line (with the possible exception of the boundary nodes, which are simply always zero in your case - no pseudonodes or anything). For example,

bx(1,t+dt)=0;bx(end,t+dt)=0;&#37;set boundary nodes
bx(2:end-1,t+dt)=dt*(i*x.*bx(2:end-1,t)+i*ux(2:end-1,t))+bx(2:end-1,t);%internal nodes

edit: This can also be written in a single line:
bx,t+dt)=[0,dt*(i*x.*bx(2:end-1,t)+i*ux(2:end-1,t))+bx(2:end-1,t),0]';
It's transposed (the apostrophe) since you're using rows for spatial nodes and columns for temporal nodes. It will also improve your speed considerably if you declare the arrays before assigning them to avoid memory fragmentation. This was a major problem when I first started programming in ML and actually makes a very large difference in performance for large problems. A simple call, such as
bx=zeros(SpatialNodes,TimeNodes)
will initialize it as a matrix of the size you need.