Solving coupled Diff Eqs with Runge Kutta

[Biftheunderstudy]

I have a set of 6 coupled differential equations, the variables are U_x, U_y, U_z, b_x,b_y,b_z and the differential is over time. For brevity I won't post the equations themselves but they are effectively:

db/dt = x*b + i*u
and
du/dt = x*u + i*b +dQ/dx

Where the bold means vector and i is the imaginary number.

Normally this is an easy application of Runge Kutta in complex variables however there is a problem.
Q is defined as an integral which has to be numerically integrated.

Q = \int (du_y/dx + i*u_x) * exp(-|x-x'|) dx'

note the absoltue value, so really its 2 integrals from 0 -> infinity and -infinity -> 0

I have implemented a finite difference method to compute the derivatives: ( u_y(2) - u_y(1) )/ dx and a similar formula for dQ/dx. All this means that my code has to carry around u_y(2), u_y(1), and u_y(0) to calculate all the derivatives for each time step.

Now where I am confused is the step where the Runga - Kutta routine is computing the various K values, heres a snippet of the code in Fortran 90:

subroutine RK4(bx,by,bz,ux,uy,uy_1,uy_2,uz)
IMPLICIT NONE
complex :: bx,by,bz,ux,uy,uy_1,uy_2,uz
complex :: k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44,k51,k52,k53,k54,k61,k62,k63,k64

k11 = .5*dt*Eq1(bx,by,bz,ux,uy,uz)
k21 = .5*dt*Eq2(bx,by,bz,ux,uy,uz)
k31 = .5*dt*Eq3(bx,by,bz,ux,uy,uz)
k41 = .5*dt*Eq4(bx,by,bz,ux,uy,uy_1,uy_2,uz)
k51 = .5*dt*Eq5(bx,by,bz,ux,uy,uy_1,uz)
k61 = .5*dt*Eq6(bx,by,bz,ux,uy,uy_1,uz)

k12 = .5*dt*Eq1(bx+k11,by+k21,bz+k31,ux+k41,uy+k51,uz+k61)
k22 = .5*dt*Eq2(bx+k11,by+k21,bz+k31,ux+k41,uy+k51,uz+k61)
k32 = .5*dt*Eq3(bx+k11,by+k21,bz+k31,ux+k41,uy+k51,uz+k61)
k42 = .5*dt*Eq4(bx+k11,by+k21,bz+k31,ux+k41,uy+k51,uy_1,uy_2,uz+k61)
k52 = .5*dt*Eq5(bx+k11,by+k21,bz+k31,ux+k41,uy+k51,uy_1,uz+k61)
k62 = .5*dt*Eq6(bx+k11,by+k21,bz+k31,ux+k41,uy+k51,uy_1,uz+k61)

k13 = dt*Eq1(bx+k12,by+k22,bz+k32,ux+k42,uy+k52,uz+k62)
k23 = dt*Eq2(bx+k12,by+k22,bz+k32,ux+k42,uy+k52,uz+k62)
k33 = dt*Eq3(bx+k12,by+k22,bz+k32,ux+k42,uy+k52,uz+k62)
k43 = dt*Eq4(bx+k12,by+k22,bz+k32,ux+k42,uy+k52,uy_1,uy_2,uz+k62)
k53 = dt*Eq5(bx+k12,by+k22,bz+k32,ux+k42,uy+k52,uy_1,uz+k62)
k63 = dt*Eq6(bx+k12,by+k22,bz+k32,ux+k42,uy+k52,uy_1,uz+k62)

k14 = dt*Eq1(bx+k13,by+k23,bz+k33,ux+k43,uy+k53,uz+k63)
k24 = dt*Eq2(bx+k13,by+k23,bz+k33,ux+k43,uy+k53,uz+k63)
k34 = dt*Eq3(bx+k13,by+k23,bz+k33,ux+k43,uy+k53,uz+k63)
k44 = dt*Eq4(bx+k13,by+k23,bz+k33,ux+k43,uy+k53,uy_1,uy_2,uz+k63)
k54 = dt*Eq5(bx+k13,by+k23,bz+k33,ux+k43,uy+k53,uy_1,uz+k63)
k64 = dt*Eq6(bx+k13,by+k23,bz+k33,ux+k43,uy+k53,uy_1,uz+k63)

bx = bx + k11/6. + k12/3. + k13/3. + k14/6.
by = by + k21/6. + k22/3. + k23/3. + k24/6.
bz = bz + k31/6. + k32/3. + k33/3. + k34/6.

ux = ux + k41/6. + k42/3. + k43/3. + k44/6.
uy = uy + k51/6. + k52/3. + k53/3. + k54/6.
uz = uz + k61/6. + k62/3. + k63/3. + k64/6.
end subroutine RK4

I feel that I should be adding K5* to uy_1 and uy_2 during this step, where uy_1 and uy_2 are uy(1) and u_y(0) for the purposes of calculating derivatives.

I can elaborate on specifics more at request.

Does anyone have any suggestions/advice?

 

[KayGee]

I have used the RK4 scheme a few times, but I really dislike having to read code. Would it be possible to see a scanned copy of the description of the problem (including the equations, boundary/initial conditions, etc.)?

 

[eLiu]

Why are you using RK4? Your error in space is O(deltax) because of the divided difference. So unless deltax is about deltat^4, RK is just wasting effort. What rule are you using to compute the integral?
Edit: and setting deltax = deltat^4 would be quite hard because first divided differences begin to fail (in double precision) around deltax = ~1e-8

As for uy_1 and uy_2: when you use RK type methods (1 step composed of several fractional steps), a good rule is to always set your BCs before you have to use them again. So before calculating k41, take uy_1 and uy_2 from the most recent time (e.g. initial data on the first step). Before calculating k42, find the new values for uy_1 and uy_2 based on the intermediate step from k*1. What you're doing now is: 1) calculate k41 by taking uy_1 and uy_2 from the most recent time, Ti 2) use those boundary values from Ti at all intermediate times. 3) update solution to time Ti+dt. 4) repeat. You cannot use the boundary values for Ti at all intermediate times; they must be recalculated with the RK method, just like all the other state variables or it's like applying Forward Euler to the BC data.

Edit2: Doing what you're doing now will -not- cause it to break. But it will likely reduce your time accuracy to 1st order since your BC data is only resolved to first order accuracy in time.

Also, F90? You're crazy Fortran makes my head hurt.

 

[CycloWizard]

I could be wrong because it's late on a Friday night and I've been drinking a little, but isn't your equation for u a PDE? If so, RK doesn't apply. If it were me, I would just use a finite difference method and be done with it. FD's are very straightforward, even for coupled PDEs, though depending on the domains they aren't necessarily very efficient.

 

[eLiu]

Originally posted by: CycloWizard
I could be wrong because it's late on a Friday night and I've been drinking a little, but isn't your equation for u a PDE? If so, RK doesn't apply. If it were me, I would just use a finite difference method and be done with it. FD's are very straightforward, even for coupled PDEs, though depending on the domains they aren't necessarily very efficient.

Huh? Suppose I want to solve a convection problem on a square domain. For completeness, let's say the velocity field swirls st the boundaries are stationary (so no BCs needed). On the interior, I have a lot of -higher order- choices available to me. If I want to consider point values FD-style, then there's say ENO/WENO, spectral, or even just higher order upwind discretizations. Now if I couple my higher order space method to a forward-euler type of time stepping, well then I just wasted a lot of effort for nothing b/c forward euler (unless I pick an obscenely small time step) will knock my accuracy back down to 1st order. Hence where RK4 comes in.

He's got the same kind of situation here. The problem is hyperbolic-ish (looks very similar to a hyperbolic problem...) with somewhat complex source terms. So the approach will be to apply some sort of space discretization to the space differential terms & then time-step.

But since he's using (u2-u1)/dx in space, his spacial accuracy is already so poor that RK4 is unnecessary. It can still be applied; it's just overkill w/o a better space discretization to go along.

Well that's my understanding anyway,
-Eric

edit: Ok now that we're here, oh god I hate FD methods. They're so nonsensical. Point values? That's totally bizarre from a mathematical point of view. I mean getting FD to work for complex domains, shock problems, etc is so hack-ish. It's quick but it's pretty damn dirty. Finite Elements is the way to rock with their mathematically rigourous variational formulations.

 

[KayGee]

Originally posted by: CycloWizard
I could be wrong because it's late on a Friday night and I've been drinking a little, but isn't your equation for u a PDE? If so, RK doesn't apply. If it were me, I would just use a finite difference method and be done with it. FD's are very straightforward, even for coupled PDEs, though depending on the domains they aren't necessarily very efficient.

It's not uncommon to use the RK4 or even RK6 scheme when solving the Euler equations associated with aeroacoustics applications. There are published papers and I personally have used it and have gotten extremely good results.

 

[Biftheunderstudy]

Ok, Thanks for all the replies.

Originally posted by: CycloWizard
I could be wrong because it's late on a Friday night and I've been drinking a little, but isn't your equation for u a PDE?

I thought that for about a week until I wrapped my head around it, but no its just a time derivative.

The other complication is that the solution should blow up at infinity so I must keep the solution truncated. The initial conditions are a flat homogeneous u and b with a gaussian bump in b.
Does anyone know how to post a pdf? I have the equations written out in LaTEX.

 

[Biftheunderstudy]

I realize that RK4 is kind of overkill and that the 1st order Finite Difference derivative is killing my accuracy, but its more of a proof of concept at this point. If it works I'll probably use a 2nd order central diffence scheme -- still losing accuracy but I can live with it. As for the integration method I'm using Simpsons rule.

 

[CycloWizard]

Originally posted by: Biftheunderstudy
I thought that for about a week until I wrapped my head around it, but no its just a time derivative.

Ah, I see now. I suppose I should skip this forum under such circumstances in the future.

The other complication is that the solution should blow up at infinity so I must keep the solution truncated. The initial conditions are a flat homogeneous u and b with a gaussian bump in b.
Does anyone know how to post a pdf? I have the equations written out in LaTEX.

You just need to find a host, or you could even paste in the TeX here.

I realize that RK4 is kind of overkill and that the 1st order Finite Difference derivative is killing my accuracy, but its more of a proof of concept at this point. If it works I'll probably use a 2nd order central diffence scheme -- still losing accuracy but I can live with it. As for the integration method I'm using Simpsons rule.

There are higher order difference schemes which use an increasing number of points. Mathworld gives some relatively confusing breakdowns. Abramowitz and Stegun give higher-order difference equations directly, but I left mine in the office, so I don't have it over the weekend.

 

[eLiu]

Originally posted by: KayGee

Originally posted by: CycloWizard
I could be wrong
snip
sufficient.

It's not uncommon to use the RK4 or even RK6 scheme when solving the Euler equations associated with aeroacoustics applications. There are published papers and I personally have used it and have gotten extremely good results.

haha, yep. The higher your space order, the better RK you'd need. Though practically there are not many people working with spacial methods of more than 5th or 6th order, esp in 2 or 3D. This is b/c it's rather difficult to form higher order discretizations using FD or FV (except for the spectral folks). You have to come up with special rules and there's a whole slew of stability issues w/problems that are interesting. But I do discontinuous galerkin finite elements, where higher order is quite natural. And like any time dependent problem (convection-diffusion, shallow water, traffic, euler, navier stokes, burger's, KdV) will require some kind of time-stepping. Even steady-state problems are typically time-stepped (start from a guess, use a coarse implicit step to converge to the steady state), b/c when you're looking at like a billion unknowns, it's impractical & inaccurate to try and solve the system directly or even iteratively.

Originally posted by: Biftheunderstudy
I realize that RK4 is kind of overkill and that the 1st order Finite Difference derivative is killing my accuracy, but its more of a proof of concept at this point. If it works I'll probably use a 2nd order central diffence scheme -- still losing accuracy but I can live with it. As for the integration method I'm using Simpsons rule.

Ok; just wanted to make sure you were aware of the issue with using discretizations of varying accuracy.

Originally posted by: Biftheunderstudy
I thought that for about a week until I wrapped my head around it, but no its just a time derivative.

Wait it's not a PDE? You have the derivative of an integral, so that doesn't change the behavior of the integrand, I thought. So the eqn is basically like u_t = f(u) + g(u)*u_x (where underscore means derivative), which is hyperbolic.

Originally posted by: CycloWizard
There are higher order difference schemes which use an increasing number of points. Mathworld gives some relatively confusing breakdowns. Abramowitz and Stegun give higher-order difference equations directly, but I left mine in the office, so I don't have it over the weekend.

It's actually extremely easy to generate FD stencils for any derivative in 1D. Then you can use tensor products to extend to higher dimensions as needed (well, assuming your higher dimensional grid is a tensor product grid!). Here's soem simple matlab code to do it... I'll let you work out why it works.

nderiv is the order of the derivative that you want. x is the points you'd like to use in the stencil. I usually specify x in multiples of h. So suppose I want to generate the classic 2-point, 1st derivative: x=[0,1], nderiv = 1 results in s -> [-1,1], meaning that your stencil will be: (u_{j+1} - u_j)/h

3 points, 1st deriv, left-biased: x=[-2,-1,0], nderiv = 1 results in s->[-2,0.5,1.5], so the approximation is: (-2*u_{j-2} + 0.5*u_{j-1} + 1.5*u_j)/h, which is more commonly written (-4u_{j-2} + u_{j-1} + 3u_j)/(2h)

Notice that you could also do like x = [-1, 0, 0.5] for a non-equidistant stencil.

V = rot90(vander(x)); % Vandermonde matrix
b = zeros(size(x))'; b(nderiv+1) = factorial(nderiv);
s = V\b; % stencil weights

Try it out & see if you can regenerate the results from the mathworld page. (vander is the vandermonde matrix. rot90 is b/c matlab's vandermonde isn't the definition I'm used to.) If you have trouble using that snippet, pm me. If you're unsure of what the error will be, you can check it like this:
s'*x.^a, where a is the power of the term you want. So for a 2nd order accurate method, you'd expect that expression to be 0 for a = nderiv+1 and nderiv+2.

If you prefer Maple, you can do:
V:= Transpose(Matrix(3,shape=Vandermonde[<-2*h,-h,0>]); #vandermonde matrix
b:=Vector([0,1,0,0]); #1st derivative
s := LinearSolve(V,b);
Then check errors: DotProduct(s,subs(a=2,Vector([(-2*h)^a,(-h)^a,0])));
would give you the error on the h^2 term.

Note: you need to be VERY careful in how you select your FD stencils! Just because two stencils are say O(h^4) DOES NOT mean that they will give the same performance. (As a classic example, FTCS--Forward Time, Centered Space--fails miserably for hyperbolic problems; but it works fine for elliptic problems.) You will have to do some kind of upwinding OR use centered space with artificial viscosity (e.g. Lax-Wendroff) for hyperbolics.

 

[Biftheunderstudy]

Bear with me, waiting to get the Tex file from my office computer.
The derivative is probably not that important in the end, so getting too complex with it is not needed. That said, I'm not very good with Matlab and only passable at Maple and I've never heard of a vandermonde
The main goal is to try and figure out how to carry around the values of u_y to compute the derivative so that RK4 works nicely.
I had initially thought to do straight FD methods but eventually changed my mind since I know RK4 better.
Anyhow, I'll likely get some more answers on Monday.

 

[eLiu]

oh a MxN vandermonde matrix will have the terms of M power series (N terms each) arrayed out in rows (or cols after you rot90 it).

FD stencils come from taylor expansions. So using vandermonde, you can build the taylor expansions for the stencil points that you're interested in. Then by choosing the right-hand-side appropriately, you do a linear solve to find what combination of stencil points (their taylor expansions really) will yield like, a 3rd derivative with no error in the 4th and 5th deriv terms (so 3rd order accurate).

 

[Biftheunderstudy]

Ok here's the pdf with the equations and a definition of Q.

Equations

oh a MxN vandermonde matrix will have the terms of M power series (N terms each) arrayed out in rows (or cols after you rot90 it). FD stencils come from taylor expansions. So using vandermonde, you can build the taylor expansions for the stencil points that you're interested in. Then by choosing the right-hand-side appropriately, you do a linear solve to find what combination of stencil points (their taylor expansions really) will yield like, a 3rd derivative with no error in the 4th and 5th deriv terms (so 3rd order accurate).

Ah, I understand. Probably overkill for my problem any way. I'll likely use a 2nd or 3rd order FD method for the derivative in the end.

So in the mean time can I use RK4 on the equation for u since it's hyperbolic? I have some experience with Godunov and Lax Wendroff but I only know how to do it for a Riemann type problems with conservative Euler equatons.

 

[racolvin]

Good lord

The fact that you guys understand all that gibberish is completely beyond me. I bow to your superior intellect

 

[eLiu]

Originally posted by: Biftheunderstudy
Ok here's the pdf with the equations and a definition of Q.

Equations

oh a MxN vandermonde matrix will have the terms of M power series (N terms each) arrayed out in rows (or cols after you rot90 it). FD stencils come from taylor expansions. So using vandermonde, you can build the taylor expansions for the stencil points that you're interested in. Then by choosing the right-hand-side appropriately, you do a linear solve to find what combination of stencil points (their taylor expansions really) will yield like, a 3rd derivative with no error in the 4th and 5th deriv terms (so 3rd order accurate).

Ah, I understand. Probably overkill for my problem any way. I'll likely use a 2nd or 3rd order FD method for the derivative in the end.

So in the mean time can I use RK4 on the equation for u since it's hyperbolic? I have some experience with Godunov and Lax Wendroff but I only know how to do it for a Riemann type problems with conservative Euler equatons.

Just out of curiosity, where does this equation come from?

As for the hyperbolic-bit, RK4 is fine. In fact the time-stepper you use for hypebolic problems practically doesn't matter. What does matter is the way you discretize the space derivative. Godonuv only makes sense in the finite volume sense; it's finite difference parallel is 1st order, -conservative- space discretizations. Lax-Wendroff has meaning in FD & FV; it's a taylor-series trick to "create" both a 2nd order time step & 2nd order space discretization at once making a correction to the classic (but awful) "forward [euler] time, centered space" scheme. The correction is derived from expansions but it can be viewed as a form of artificial viscosity.

Anyway, back to your problem: after going through your derivation, I think your problem is one of reverse-spreading? Like backwards diffusion but it won't go like exp(-x^2), it'll go like exp(-x). Maybe? I'm not completely sure b/c my work doesn't involve sqrt(-1) in my governing equations, aahhh! If that's the case, then conservative vs non-conservative discretization isn't a big deal b/c you don't have shocks (discontinuities in the soln), is that right? Then with the factors of i on all your space derivatives, I think you will not have to worry about upwinding either (but i could be wrong about that). (That is, advection equation has the dispersion relation w ~ ik, so in the fourier domain that's just a time-shift... but your eqn will have w ~ -k, which is exponential growth or decay.) If that's the case, then you can do anything for your space discretization. Once you move past 1st order upwind, I'd stick with stencils that are centered in space. (Again this is all assuming your problem isn't advective; if so, you will want to avoid centered/symmetric stencils & use upwinding/downwinding as needed.)

Stop me if I'm getting something wrong here...

-Eric

racolvin, eh for me, it's just been a lot of classes. So if you want to take some courses in numerical analysis, linear algebra, and PDEs, you too could be part of this party! haha but possibly you have more exciting things to do

 

[Biftheunderstudy]

The equations are from MHD.
I'll be moving to centered space derivatives to keep things as symmetrical as possible, but I don't know about going to higher order. The derivative means I have to carry 3 uy's around which means that my RK4 is already gigantic. So I have to do 3 steps of RK4 at once which without subroutines is 72 lines of intermediate steps :S
Going to higher order is going to make that pretty crazy. For now first order in space I think I can live with.

 

[eLiu]

MHD? Google leads me to... magnetohydrodynamics? sounds tough. (also maybe the MA highway dept)

So if you're moving to centered, I take it that what I said in the post above is more or less correct? In that case, my calling the problem hyperbolic at first was an error on my part... oops. The sqrt(-1) factors threw me off until I tried to apply a little fourier analysis.

Ok well good luck with your project, and good luck using Fortran, blech

 

[Biftheunderstudy]

Yup Magnetohydrodynamics is right...also tough is right. Calling them hyperbolic is not far off the mark since MHD is VERY closely related to fluid dynamics, there are just some other things to worry about. In my case, the fluid is incompressible so I don't have to deal with shock solutions...yet. Your comments in your post were eerily echoed by my supervisor.
In any case, I now have to wrap my head around how to make subroutines and functions so that I can cut down the number of lines I need to make the code at least partially readable. Fortran 90 in this case is actually a godsend with the array handling, I cringe at the thought of trying to do it in C. Fortran is still alive and kicking in the astrophysics community, though I would like to sit down and learn some more C++ so I can have something useful when I graduate.

 

[eLiu]

No shocks is a plus. Truth be told shock capturing with finite differences is somewhat painful. As for the hyperbolicity, formally hyperbolic problems cannot have complex eigenvalues when you write them (or their linearized forms) in u_t + A*u_x = RHS form. So technically I think your system isn't hyperbolic; the i factors turn traveling waves into growth/decay modes, I think.

I code in C... a Fortran guy telling me he's got the better deal? Now I've seen it all, haha Yeah I know F is alive and kicking. I rarely code in it, but I interface with the BLAS and LAPACK all the time. There are still plenty of people/codes that are mostly/completely F. But I'd agree that it couldn't hurt to pick up some C or C++.