52 Reputation

4 Badges

12 years, 301 days

MaplePrimes Activity

These are answers submitted by Benj

This actually turns out to be pretty workable.  The 'op' command is the key that I was missing, the maple help never really helped there.  evalf cleans it up a bit, this whole thing is supposed to be an approximation to a larger problem so I'm just looking for a mostly approximate answer.

My next question though, is that the initial conditions don't seem to hav worked that well.  (I'm guessing it is because omega and x are paramaters that haven't been defined yet)

The reason I say this, is that when I evaluate the solution at t=0 I would expect to get my initial conditions right?  If not is there an explanation for this?

Thanks again.

I'm sorry I made a mistake in my posting, the correct equations should be:

 eq1 := diff(RBx(x, t), t) = -x*IBx(x, t)-IUx(x, t)
eq2 := diff(IBx(x, t), t) = x*RBx(x, t)+RUx(x, t)
eq3 := diff(RBy(x, t), t) = -x*IBy(x, t)-IUy(x, t)-3*RBx(x, t)/(2*omega)
eq4 := diff(IBy(x, t), t) = x*RBy(x, t)+RUy(x, t)-3*IBx(x, t)/(2*omega)
eq5 := diff(RBz(x, t), t) = -x*IBz(x, t)-IUz(x, t)eq6 := diff(IBz(x, t), t) = x*RBz(x, t)+RUz(x, t)
eq7 := diff(IUx(x, t), t) = -x*IUx(x, t)-IBx(x, t)+2*RUy(x, t)/omega
eq8 := diff(RUx(x, t), t) = x*RUx(x, t)+RBx(x, t)+2*IUy(x, t)/omega
eq9 := diff(RUy(x, t), t) = -x*IUy(x, t)-IBy(x, t)-RUx(x, t)/(2*omega)
eq10 := diff(IUy(x, t), t) = x*RUy(x, t)+RBy(x, t)-IUx(x, t)/(2*omega)
eq11 := diff(RUz(x, t), t) = -x*IUz(x, t)-IBz(x, t)
eq12 := diff(IUz(x, t), t) = x*RUz(x, t)+RBz(x, t)


> sys := [diff(Bx(t), t) = I * x*Bx(t) + I * Ux(t),
          diff(By(t), t) = I * x * By(t) + I * Uy(t) - 3 * Bx(t)/(2*omega),
          diff(Bz(t), t) = I * x * Bz(t) + I * Uz(t),
          diff(Ux(t), t) = I * x * Ux(t)  + I * Bx(t) - 2 * I * Uy(t)/omega, 
          diff(Uy(t), t) = I * x * Uy(t) + I * By(t) - Ux(t)/(2*omega),
          diff(Uz(t), t) = I * x * Uz(t) + I * Bz(t)];

I've made it to this step but the eigenvectors are a huge mess of RootOf things and wierd constants.  Is there an easy way to get this to a nice looking solution including the initial conditions?

If I try to just pdsolve(sys) it gives me a bit nicer solution but I can't figure out how to set the initial conditions in pdsolve.  If I define a list with all the initial conditions on it like this:

IC := [Bx(0) = exp(-x^2/(2*sigma^2)), By(0) = 0, Bz(0) = (3*I)*k*(diff(exp(-x^2/(2*sigma^2)), x))/(2*omega), Ux(0) = 0, Uy(0) = 0, Uz(0) = 0]

and then call pdsolve(sys,ICs), I get this error:

Error, (in pdsolve/sys) too many arguments; some or all of the following are wrong: [{Bx(t), By(t), Bz(t), Ux(t), Uy(t), Uz(t)}, [Bx(0) = exp(-(1/100)*x^2*ln(2)), By(0) = 0, Bz(0) = -((3/100)*I)*k*x*ln(2)*exp(-(1/100)*x^2*ln(2))/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0]]

How might I fix this?

Here are the fully expanded equations written in LaTex, should be able just to copy and paste into a Tex editor to decipher.  The dependance of these equations should all be on x so hopefully Maple won't have too much difficulty.

After a bit of manipulation, the real and imaginary parts can be separated yielding:

\[\partial_t b_{xr} = -x b_{xi} - u_{xi}\]
\[\partial_t b_{xi} = x b_{xr} + u_{xr}\]

\[\partial_t b_{yr} = -x b_{yi} -\frac{3}{2 \bar{\omega}_A} b_{xr} -u_{yi}\]
\[\partial_t b_{yi} = x b_{yr} -\frac{3}{2 \bar{\omega}_A} b_{xi} +u_{yr}\]

\[\partial_t b_{zr} = -x b_{zi} -u_{zi} \]
\[\partial_t b_{zi} = x b_{zr} +u_{zr} \]

\[ \partial_t u_{xr} = - x u_{xi} +\frac{2}{\bar{\omega}_A} u_{yr} - \frac{3}{2 \bar{\omega}_A} \partial_x Q_r - b_{xi} \]
\[ \partial_t u_{xi} =  x u_{xr} +\frac{2}{\bar{\omega}_A} u_{yi} - \frac{3}{2 \bar{\omega}_A} \partial_x Q_i + b_{xr}\]

\[ \partial_t u_{yr} = - x u_{yi} -\frac{1}{2 \bar{\omega}_A} u_{xr} - Re(iQ) - b_{yi} \]
\[ \partial_t u_{yi} =  x u_{yr} -\frac{1}{2 \bar{\omega}_A} u_{xi} - Im(iQ) + b_{yr} \]

\[ \partial_t u_{zr} = - x u_{zi} - Re(\frac{i}{\kappa} Q) - b_{zi}\]
\[ \partial_t u_{zr} =  x u_{zr} - Im(\frac{i}{\kappa} Q) + b_{zr}\]

Where the subscript $r$ and $i$ refer to the real and imaginary parts respectively.
The real part of $\partial_t Q$ is given by:
\[ \frac{4}{3} u_{yr} - \int_x^{\infty} \left( \frac{2}{3} k u_{yr} - \frac{2}{9} \bar{\omega}_A u_{xi} \right) e^{k(x-x')} dx' - \int_{-\infty}^x \left( \frac{2}{3} k u_{yr} + \frac{2}{9} \bar{\omega}_A u_{xi} \right) e^{-k(x-x')} dx'\]
And the imaginary part is given by:
\[ \frac{4}{3} u_{yi} - \int_x^{\infty} \left( \frac{2}{3} k u_{yi} + \frac{2}{9} \bar{\omega}_A u_{xr} \right) e^{k(x-x')} dx' - \int_{-\infty}^x \left( \frac{2}{3} k u_{yi} - \frac{2}{9} \bar{\omega}_A u_{xr} \right) e^{-k(x-x')} dx'\]
The real part of $iQ$ is given by:
\[ \int_x^{\infty} \left( \frac{2}{3} u_{yi} + \frac{2}{9 k} \bar{\omega}_A u_{xr} \right) e^{k(x-x')} dx' + \int_{-\infty}^x \left( -\frac{2}{3} u_{yi} + \frac{2}{9 k} \bar{\omega}_A u_{xr} \right) e^{-k(x-x')} dx' \]
And the imaginary part of $iQ$ is:
\[ \int_x^{\infty} \left( -\frac{2}{3} u_{yr} + \frac{2}{9 k} \bar{\omega}_A u_{xi} \right) e^{k(x-x')} dx' + \int_{-\infty}^x \left( \frac{2}{3} u_{yr} + \frac{2}{9 k} \bar{\omega}_A u_{xi} \right) e^{-k(x-x')} dx'\]

Sure, here are the equations in expanded form

\partial_t b_x = i x b_x + i u_x
\partial_t b_y = i x u_y -\frac{3}{2 \bar{\omega}_A} b_x + i u_y
\partial_t b_z = i x b_z + i u_z
\partial_t u_x = i x u_x + \frac{2}{\bar{\omega}_A} u_y -\frac{3}{2 \bar{\omega}_A} \partial_x Q + i b_x
\partial_t u_y = i x u_y - \frac{1}{2 \bar{\omega}_A} u_x - i Q + i b_y
\partial_t u_z = i x u_z - \frac{i}{\kappa} Q + i b_z

I'm running a blank as to how I should split up these equations into real and imaginary parts given that b and u are both complex.

Also, I forgot to mention before that the variable k is different from \kappa.  k is used to solve the helmholtz equation and is given by:

k = \frac{2}{3} \bar{\omega}_A \sqrt{ 1 + \frac{1}{\kappa^2} }

I have already split the equations into components, but for ease of writing them here I left them in a more compact form.  As for splitting into real and imaginary parts, I left them as is since Fortran does complex numbers and it was easier than making my RK4 routine twice as big.  But in principle I can do that.

My main concern however is setting this problem up in Maple, primarily the syntax and dependancies that I'll need.

Page 1 of 1