Differential Equations with constraints - Lagrange multipliers

[Biftheunderstudy]

Once again digging up old math problems. I have a set of equations that look roughly like the incompressible Navier-Stokes equations. (Wish this forum supported Tex...)

\partial_t u = RHS + \nabla P

\nabla \cdot u = 0

I haven't actually tested this in LaTex so don't know if it works, here it is in words...

Its a set of 1st order linear equations with a gradient of a pressure on the right hand side. The second equation is simply the divergence of the u-field is equal to zero. The pressure acts as a lagrange multiplier here to enforce the divergence free condition.

Currently, I've been using an implementation of the projection method used for the NS equations but I'm worried about the accuracy since I'm using a 1st order method to solve the Poisson equation.

What I'm looking for is another way to solve this type of system. I have access to Maple/Matlab/Mathematica etc. ... just not the knowledge to solve this type of system using them.

Those equations are of course a simplification, u is a vector and there are 3 other coupled equations to go with them so I would need a vectorized solution method as well.

 

[CycloWizard]

That should work in TeX, and I happen to be solving a very similar system in MATLAB at this very moment (browsing ATHT while it solves...). MATLAB can solve this type of problem pretty easily depending on the dimensionality of the problem, among other things. What is "RHS"? I assume it's some sort of Laplacian.

Assuming you have only one spatial dimension, you can use the pdepe function to solve the PDE. This is actually easier to implement than the ODE solvers for boundary-value problems (at least for me). You essentially compute three parameters c, f, and s (corresponding to the form c \partial_t u= \partial_x (f)+s, so f is the flux term and s is the source term). The output from your PDE function is then just [c,f,s] as column vectors, where each row is the parameter for one equation in your system. Let me know if you have any questions or want some code/examples. I'm not sure if this one can handle solving for parameters like Lagrange multipliers or not, as I haven't had to worry about that.

For the steady solution, I'd recommend reading this for an overview, then the documentation for the bvp4c solver (though I should note that there is also now a bvp5c solver in newer versions). I know this one can solve for parameters and that one of the bvp4c documentation examples gives just such a case. The notation for these solvers is a bit tricker, as you have to recast the problem as a system of first-order ODEs, but I can help with this if you need it (it's easy once you play around with it a bit).

 

[Born2bwire]

I don't do work with CFD but I do computational electromagnetics and most of the techniques are the same. It seems like you could easily do a finite difference time domain algorithm for the problem and use Matlab to step through the solution. The codes are fairly simple to write but it's hard to say if that is applicable as opposed to solving the differential equations directly without knowing your boundary conditions. The FDTD solution is nice for the fact that you can have complicated boundary conditions and it doesn't affect the complexity of the coding or the runtime of the problem. Though of course it would be overkill for simple problems like a cavity.

 

[eLiu]

We need some more information... most importantly, what shape is the domain? How accurate do you need to be? How big (spatial & temporal degrees of freedom) is the problem? Boundary conditions would also be nice but it's not really that relevant.

Also what are the other 3 equations? I see 4 equations here with 4 unknowns (3 components of "velocity" and "pressure").

And where are you needing to solve a poisson problem...? The equation you wrote is not elliptic. It's more like the incompressible euler equations with the (nonlinear) convective term removed. (Well, depends on what "RHS" looks like.)

And lastly, how the heck did someone manage to write a *first order* poisson solver? Even the most naive finite difference implementation is 2nd order accurate...

 

[Biftheunderstudy]

Hmm, ok. Tried to keep them simple but here they are:

\partial_t b_x = i x b_x + i u_x
\partial_t b_y = i x b_y + i u_y - 3/(2 w) b_x
\partial_t b_z = i x b_z + i u_z
\partial_t u_x = i x u_x + i b_x -2/w u_y - 3/(2 w) \partial_x Q
\partial_t u_y = i x u_y + i b_y - 1/(2 w) u_x - i Q
\partial_t u_z = i x u_z + i b_z - i/\kappa Q
3/2 \partial_x u_x + i w u_y + i w/kappa u_z = 0

Those are the equations, the initial conditions are:
b_x(x,0) = e^{-x^2 / (2 \sigma^2)}
b_y(x,0) = 3 i / (2 w) \partial_x b(x,0)
b_z(x,0) = u_x(x,0) = u_y(x,0) = u_z(x,0) = 0

I am trying to solve this on a domain from either -50 to 50 or 0 to 100. The latter case requires altering the gaussian to peak at x=50. sigma is defined such that the full width half maximum of the gaussian is around 10.

The subscript _x, _y, _z denote the component of that vector i.e. b_x is the x component of the vector b.

\partial_t and \partial_x are partial derivatives in time and space respectively.

The last equation is the divergence equation which acts as a constraint and is where the difficulty comes in. Q is like a pressure and must be solved such that the divergence of u stays zero. (Note that the equations are in a form such that the nabla operator is (1.5*\partial_x, i*w, i*w/kappa) ).

I've done a finite difference approach to this using a projection method (Matlab runs out of memory before the solution becomes stable due to keeping the arrays of state variables around). There is a good page on the projection method on wiki...

In the end I'm using an RK4 method to solve the 6 equations without the pressure term (Q) and then correcting it by solving Q" - K^2 Q = 9/4w * div u / dt. I'm using a backwards euler method to solve this since I can't get it to solve stably using an explicit method. Given the circumstances I can't get a higher order implicit solver to work properly and so I've been stuck with a low order Euler approach.

I've tried a few things with pdepe but I'm not sure how to go about it for this set up. Can pdepe solve the incompressible momentum equation? If so, it can probably solve this.

 

[CycloWizard]

Ah, now I see why you didn't write out RHS originally. I believe pdepe should be able to handle this scenario as follows:

%Main script:
NumXNodes=100;%number of spatial nodes
NumTNodes=100;%number of temporal nodes
m=0;%coordinate symmetry (0=Cartesian geometry)

%Generate the mesh
x=linspace(-50,50,NumXNodes);%generate the mesh
t=linspace(0,t_max,NumTNodes);%temporal nodes (may use logspace here instead for more efficient solution - usually helps in my problems)

%Solve the problem
sol = pdepe(m,@pdefun,@pdeinit,@pdebc,x,t);

%Extract the solution values
b_x=sol(end,:,1);
b_y=sol(end,:,2);
%and so on for other variables

%Functions
function [c,f,s]=pdefun(x,t,u,DuDx)%defines the PDE parameters

c=ones(7,1);%coefficients of the time derivative
f=[0,0,0,Q,0,0,3/2*u_x];%"flux" term (anything acted on by d/dx)
s=[i*x*b_x+i*u_x...et cetera (all of the leftover terms in the other equations which I'm not going to bother typing out)

function [b]=pdeic(x)%initial conditions
sigma=20;%whatever value it needs to be for FWHM=10?
b=[e^(-x^2 / (2 \sigma^2))%b_x
	3 i / (2 w) \partial_x b(x,0)%b_y
	zeros(4,1)];%b_z=u_x=u_y=u_z=0
	
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)%boundary conditions
	%specify your BCs - these are all Dirichlet (i.e. b=u=0 at both boundaries)
	pl=ul;%left-most boundary (x=-50)
	pr=ur;%right-most boundary (x=+50)
	
	%Specify no Neumann BC contributions
	ql=zeros(size(ul));%x=-50
	qr=q;%x=+50

edit: Don't have MATLAB on this computer, so I can't test this out for you. Anyway, it should be enough to get you started and keep you out of trouble for a little while at least.

 

[Biftheunderstudy]

Ok, so first question. Won't this matlab implementation simply pick out the trivial Q = 0 solution and then solve \partial_x u_x= -i w u_y -i w/kappa u_z for b_x such that div u remains 0?

Next, my matlab is very rusty as I just picked it up a couple months ago and I'm entirely self taught in this regard.

In:
c=ones(7,1) for the time derivatives shouldn't this be
c = ones(6,1) ;
c(7) = 0; %the last equation has no time derivatives?

Next,
f=[0,0,0,3/(2*w)*Q,0,0,3/2*u_x];
u needs to be sent to pdefun as a vector so I should probably make u_x a dummy variable which is a component of the 7 part vector u right? Same with Q? i.e. u(1,: ) = b_x, u(4,: )=u_x, u(7,: ) = Q or would u(1)=b_x be more appropriate?

Lastly, the very last line:
qr=q;%x=+50

should that be qr=ql?

edit:
Also, looking at the help for pdepe, shouldn't c,f,s be column vectors? like f = [0;0;0;3/(2*w)*Q;0;0;3/2*u_x]; etc?

 

[CycloWizard]

If Q=\vec{0} is a valid solution, then MATLAB probably would find it as such. However, I think your continuity equation (or whatever you want to call your seventh equation) will ensure that Q has the correct value, as this makes it a problem of seven variables and seven unknowns if I'm understanding you correctly.

Oops, you're right about the c declaration. It will be a little faster if you do it like this:

c=ones(7,1);
c(7)=0;

MATLAB can slow way down if you don't preallocate the array at its final size, which is what I started doing there. I just forgot to change that last value.

The order doesn't matter as long as you're consistent. You can send u as a vector by setting the options for pdepe as follows:

BVPOptions=bvpset('Vectorize','on');

I'm not 100% sure whether this will work since the documentation is a bit ambiguous. It says in the bvpset help file that the 'Vectorize' command will not work with stiff solvers, and the time integration is performed using ode15s, which is a stiff solver. However, it may work for solving the spatial ODE. The vectorization you stated (i.e. u(1,) should work if the Vectorize command holds. Otherwise, I'm not really sure what will happen.

In the BCs, it should indeed be qr=ql. And yes, everything needs to be column vectors. Unfortunately, I'm also self-taught in MATLAB (indeed, I never got any training in programming at all ), so I always stumble through things unless I can test them for myself.

 

[eLiu]

godamnit, the forum ate my post. well here's a much shortened version:

I'm not sure that pdepe will work. The documentation says:
Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D

OP claims his problem is hyperbolic, not elliptic/parabolic (the incompressible euler equations are most definitely hyperbolic). The fundamental nature of these two classes (hyper vs ellip) means that *entirely* different solution methods are necessary for each. I have a good guess as to how pdepe works. They tell you the form of the PDE that it can solve in the docs. So they'll use a (symmetric) finite difference stencil to approximate the spatial derivatives in terms of the unknowns at a given time-step. This is results in a system of equations. Fourier analysis on linear elliptic problems tells us that (\delta t) has to scale like (\delta x)^n, where n is the order of the highest x-derivative, so explicit time stepping is impractical. Hence why pdepe uses ode15s (an implicit, multistep ode solver) to perform the time stepping. Solving elliptic problems goes the same way, except you don't do the time stepping bits.

In elliptic problems, if I perturb the problem on one side of the domain, the effects are felt everywhere (but in a "smoothed" out sense). Hyperbolic problems do not have this property. Instead, information travels in 1 direction only (well, 1 direction per equation). It's like wind blowing. You feel the effects of what is happening upwind of you: like if I fart upwind of you, you smell it; if I fart downwind, you'd never be the wiser. So even the most naive finite difference discretizations of hyperbolic problems must take this into account. To approximate \partial_x u, I cannot use a symmetric stencil b/c information about u only travels in 1 direction. I have to pick either u_{n} - u_{n-1} or u_{n+1}-u_{n}, depending on the direction of information travel. This is called "upwinding." pdepe will not perform upwinding.

[Small caveat: if you have a problem with a hyperbolic & elliptic parts (e.g., convection-diffusion equation), pdepe MIGHT work if the diffusive part is strong enough everywhere. There's very solid theory behind this: look up coercivity with respect to variational (finite element) methods.]

That said, I am somewhat confused about the hyperbolicity of this problem. Not going to lie, I've never even though about the (simplified) problem \partial_t u = i*u before. Fourier analysis does seem to indicate that this will have traveling wave solutions... with dispersive effects, I think? (Well "backwards" dispersion where higher wave numbers travel more slowly.) Is that accurate?

Also, sorry but I've never heard of this 'projection method' before. My training is in (discontinuous galerkin) finite elements & finite volume methods for CFD applications. (So I never have sqrt(-1) in my problems either!) At a glance, projection looks a lot like operator splitting, which is typically 1st order accurate unless you're very careful. (The fact that you're splitting the problem up kills the accuracy, regardless of how accurately you are solving the individual parts.) I don't have time now to read up on projection, maybe later.

edit: OP, do these crazy equations of yours have a name?

 

[Biftheunderstudy]

They are dimensionless, linearized MHD equations. That said the problem is fully linear and there is even an analytic solution for when there is no Q term (which I can recover).

The projection method is similar to operator splitting, its formal name is Helmholtz-Hodge decomposition. I couldn't find any information on accuracy but I get the impression that the method is an exact way to decouple the 2 parts of the problem

I have no idea how to tell if my set of equations are hyperbolic or parabolic. Judging that the most familiar looking relative to these equations is hyperbolic that would be my guess. It also makes sense since the problem involves a speed with wave solutions etc.

That said, they are barely even pde's in the first place hence why the projection method is so attractive since it breaks the problems into separate ode's to solve.

pdepe is doing something, I don't know that what its doing is correct yet...
I can't seem to get the resolution to anything close to accurate without running out of memory though.

Any other guesses on how to solve it? I was hoping that Maple might be able to solve it as a DAE, but again no luck there.

 

[Biftheunderstudy]

Here is the piece of matlab code that I've been using, it yells at me saying "DAE appears to be of index greater than 1"

Not sure why its giving me that...

function pdediv
global w;
global kappa;
w = 0.01;
kappa = 0.0001;
NumXNodes=2000;%number of spatial nodes
NumTNodes=100;%number of temporal nodes
m=0;%coordinate symmetry (0=Cartesian geometry)
t_max = 2.0;
%Generate the mesh
x=linspace(0,100,NumXNodes);%generate the mesh
t=linspace(0,t_max,NumTNodes);%temporal nodes (may use logspace here instead for more efficient solution - usually helps in my problems)

%Solve the problem
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);

%Extract the solution values
b_x=sol(:,:,1);
b_y=sol(:,:,2);
%b_z=sol(:,:,3);
%u_x=sol(:,:,4);
%u_y=sol(:,:,5);
%u_z=sol(:,:,6);
Qo = sol(end,:,7);
%and so on for other variables

figure
surf(x,t,log(abs(b_x)))

figure
plot(x,Qo)
%Functions
function [c,f,s]=pdefun(x,t,u,DuDx)%defines the PDE parameters
global w;
global kappa;
bx = u(1);
by=u(2);
bz=u(3);
ux=u(4);
uy=u(5);
uz=u(6);
Q=u(7);

c=ones(7,1);%coefficients of the time derivative
c(7)=0;
f=zeros(7,1);
f(4) = -3/2/w*Q;
f(7) = 3/2*ux;%"flux" term (anything acted on by d/dx)
s = zeros(7,1);
s(1)=1i*x*bx+1i*ux;
s(2)=1i*x*by + 1i*uy - 1.5*bx/w;
s(3)=1i*x*bz + 1i*uz;
s(4)=1i*x*ux + 1i*bx -2/w*uy;
s(5)=1i*x*uy + 1i*by + 0.5*ux/w - 1i*Q;
s(6)=1i*x*uz + 1i*bz - 1i*Q/kappa;
s(7)=1i*w*uy + 1i*w/kappa*uz;

function [b]=pdeic(x)%initial conditions
global kappa;
global w;
sigma=10/2/sqrt(2*log(2));%whatever value it needs to be for FWHM=10?
b=[exp(-(x-50).^2 / (2*sigma^2));...%b_x
	3i/(2*w*sigma^2)*(x-50).*exp(-(x-50).^2 / (2*sigma^2));...%b_y
	zeros(5,1)];%b_z=u_x=u_y=u_z=0


function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)%boundary conditions
	%specify your BCs - these are all Dirichlet (i.e. b=u=0 at both boundaries)
	pl=ul;%left-most boundary (x=-50)
	pr=ur;%right-most boundary (x=+50)
	
	%Specify no Neumann BC contributions
	ql=zeros(size(ul));%x=-50
	qr=ql;%x=+50

 

[CycloWizard]

eLiu may be right that pdepe will not solve the problem... As you said, it's a bit difficult to classify this system (at least for me ). You might check by trying to solve it and matching the solution to your known analytical solution where Q=0, as this should be a simple adaptation of your code. MATLAB is the only language I'll be any help with, so I hope I didn't push you in the wrong direction with it...

 

[eLiu]

They are dimensionless, linearized MHD equations. That said the problem is fully linear and there is even an analytic solution for when there is no Q term (which I can recover).

The projection method is similar to operator splitting, its formal name is Helmholtz-Hodge decomposition. I couldn't find any information on accuracy but I get the impression that the method is an exact way to decouple the 2 parts of the problem

I have no idea how to tell if my set of equations are hyperbolic or parabolic. Judging that the most familiar looking relative to these equations is hyperbolic that would be my guess. It also makes sense since the problem involves a speed with wave solutions etc.

That said, they are barely even pde's in the first place hence why the projection method is so attractive since it breaks the problems into separate ode's to solve.

pdepe is doing something, I don't know that what its doing is correct yet...
I can't seem to get the resolution to anything close to accurate without running out of memory though.

Any other guesses on how to solve it? I was hoping that Maple might be able to solve it as a DAE, but again no luck there.

Yeah I know it's linear MHD... good luck with that. I've never worked w/those equations before.

Still haven't read up much on that projection thing (just wiki). Operator splitting (in general) is not exact, but this may be special. Also (unless I'm misunderstanding your statement), any PDE discretization effectively breaks the PDE into a system of ODEs.

And in general, MATLAB is not a good choice for large problems... for things that boil down to (sparse) systems of equations, large problems is somewhere in the neighborhood of a few hundred thousand unknowns. (More if there's no time stepping, less if there is time-stepping)

Anyway I'm fairly sure your equations are hyperbolic. Like I said if you write u_t = i*u (i*x*u doesn't change anything; _x means partial deriv wrt x), unless I did it wrong, fourier analysis gives traveling wave solutions (no decay/growth) with dispersion. That said, except for very simple equations, I don't have a series of steps you can apply to classify a set of equations. What I do know is that elliptic requires a second order spatial operator (actually I think any positive even order is elliptic); parabolic is elliptic + time derivative; hyperbolic is time deriv + 1st order spatial deriv (maybe any odd order is hyperbolic?). Many equations don't strictly fit into one of these categories since you can easily mix them up (e.g., navier stokes eqns). The previous statement only applies if i (sqrt(-1)) is not present... I think what you really need to examine is what the fourier solution of the (linearized) problem looks like. I think spatial dependencies on (wave number)^(odd power) are hyperbolic and (wave number)^(even power) are elliptic? Maybe?

I'm not aware of an automatic solver for hyperbolic systems of equations in maple/matlab/mathematica. Not to mention that the symbolic languages will be slow for this task & even more ram-limited than matlab. That said, I would be surprised if you cannot find prebuilt codes online for solving these equations. There are many people working on computational MHD problems... you ought to be able to find like a C source you can download & compile.

I don't really have any implementation suggestions b/c all of my experience comes from solving conservation law type problems. i.e., problems that can be written like u_t + f(u)_x = source, where f(u) is the flux. This problem doesn't quite fit into that framework...

 

[Biftheunderstudy]

Keep in mind, _x meant the x component of the vector not the spatial derivative.

To my knowledge, no one has solved these equations before (which is kind of the point).

If the system stays at least block triagonal I can invert the matrix using sparse techniques in a real programming language much faster than Matlab. Currently I'm timestepping with RK4 and doing the pressure correction step using a backwards euler matrix inversion. I tried this in Matlab and could only push the resolution to about a factor of 10 bigger than my could could easily do. The code was many many times faster than Matlab even at 10x lower(higher?) resolution. This probably stems from the fact that I'm terrible at Matlab and had it poorly optimized but it shouldn't result in the huge performance differences seen. My only conclusion would be that Matlab is far slower than traditional programming languages as soon as you take away Matlab's ability to vectorize.

That said I am programming in Fortran 90 and its array handling is superb as well as the ability to also vectorize array operations much like Matlab does. (Plus built in complex number handling is key)

P.S.
Thats always something thats bugged me, saying higher resolution seems wrong ?!?!
i.e. .001 is higher resolution than .01

 

[Biftheunderstudy]

So if Matlab is not appropriate and going along the same vein as the projection method, does anyone have a good algorithm to solve this equation:

Q" - K^2 Q = F

where F = 4*w*w/9 * div u

and the divergence operator is defined as:

( 3/2 \partial_x, i w, i w / kappa) and i = sqrt(-1)

and K=2/3*w*sqrt(1+1/kappa^2)

u in the above equation is an intermediate result of solving the initial equations with no Q term.

Currently I'm solving this equation using a backwards euler method which is:

y(n+1) = y(n) + h*F(n+1)

In more detail this is:

[Q(n+1) - 2 Q(n) + Q(n-1)] / h^2 = K^2 Q(n) + F

By treating each node as a separate variable so at each node there is a algebraic equation you end up with a tridiagonal matrix, A, with the -1 and +1 diagonal having a value of 1 and the main diagonal having a value of -(2+ h^2 K^2).

Now, you have a matrix equation: A*Q = F*h^2 which can be solved via the Thomas algorithm.

Keep in mind that due to time stepping the original equations I have F as an array of numbers and it is not defined as a function so fractional step solvers will not work in this case. (i.e. RK4)

Also, I'm using backwards euler since forward euler proved to be unstable regardless of the stepsize I used so I suspect a implicit solver is necessary here.

 

[eLiu]

Few notes:
--Yeah sorry, I didn't mean to type u_x. u_t = i*u_x wouldn't be hyperbolic anyway (I think) since that solution admits exponential growth/decay (depending on the sign).

--Fully vectorized matlab code is typically something like 5x slower than well written C. And 'fully vectorized' entails you're actually doing something... like doing
A = rand(1000); b=rand(1000,1);
A\b
won't be much slower than doing the equivalent operations in C and calling LAPACK for the linear solve. Vectorizing is the ONLY way to make MATLAB code 'fast'. MATLAB is an interpreted language, so when you write a for loop, the interpreter is being called constantly to 'decode' the loop contents & make the appropriate compiled function calls. Vectorization allows the interpreter to make 1 function call to some compiled code (which has the loop built-in); this is order(s) of magnitude faster. If you want to use matlab, spend time learning how to vectorize properly.

--Solving tridiagonal systems is trivial. Additionally, LAPACK's tridiagonal code will be no worse than yours; no reason to re-invent the wheel. (Be sure to install some optimized BLAS with it, as opposed to the reference BLAS, if you care about speed.)

--I don't think this is the first time anyone has tried to solve the MHD equations... like this:
http://www.springerlink.com/content/m734l7h364072u81/
seems to be describing a finite-volume MHD solver & includes a survey of some existing methods. And unless there are many types of MHD equations (again, i know nothing about this field), if someone is solving the full-on problem, someone else has already worked on the linear form. (Because you need to know how to solve the linearized problem to solve the nonlinear problem.)

--RK4 WILL work for that problem. But you need to break it into a system of 2 first order equations first. Any ODE solver can solve that problem. Some may require tiny tiny time steps to remain stable due to stiffness (an explicit method will work, but it might not be practical... if it's failing, you're doing it wrong). If you think the ODE solver is the limiting factor here, then you can throw in higher order time steppers... like implicit RK4, a higher order backwards "euler" type scheme, some kind of implicit multistep method, etc. Lots of choices exist, but I doubt this is the limiting factor.

Also, in your Q'' +blah = F, is F known? Like F is 'data' and not a function of Q? Is this an IVP or a BVP? (i.e., where do you have boundary/initial conditions for Q?) If it's a BVP, then backwards euler is entirely the wrong approach anyway. Sounds like it is a BVP b/c you described it as a poisson problem. You do not use ODE solvers for initial value problems to solve boundary value problems.

 

[Biftheunderstudy]

I've been hesitant to use LAPACK since I'm not sure how it works with complex numbers and I'm not particularly fond of making my problem twice as big to break it into complex and real solutions...

While many people work on MHD equations both linear and non-linear, the particular version that I'm trying to solve hasn't been solved yet.

If you follow the projection method closely it gets to a line which is like:
-take the divergence of both sides (this is how the divergence constraint is enforced)
When you do that you normally get:
div u/dt = laplacian Q
Now simply applying the fact that the nabla operator is (3/2 \partial_x, i w , i w / kappa) it becomes a little different:
4 w/ 9 div u/dt = Q" - K^2 Q
Thus to apply the correction term I need to solve this.

I guess it probably could be a BVP with Q going to zero at the left and right boundaries. In fact while I was looking for a method to solve this I used bvp4c from Matlab and it matched the result from my integration and backward euler. In the end I didn't know how to write my own bvp solver that was guaranteed to work and it was much to slow in Matlab.

Now, the div u part is calculated from u which is stored in an array from the time stepping and would play the role of a source term. Now hopefully you can see why I said Laplacian in the above post, but when I apply it in my case its more like a Helmholtz equation with a source term.

Why I said RK4 won't work right out of the box is that there is a step like:
K2 = h*f(x+h/2,y+K1/2)
and my source term (f) is defined only integer mesh points while the RK4 method is asking for a value at x+h/2 (in between the grid points). Now I know I can interpolate but it adds a lot of overhead on the problem and I don't know if thats the right way to go about it.
As for the stiffness, the resolution is pretty tiny already but not enough to get forward euler working.

It also has a solution given by an integral which could be solved as well though it takes way longer than solving the ode for some reason.

Now for the time being, you can assume my RK4/simpson whatever implementation is correct, I've double/triple/quadruple checked against test problems to ensure they are working as intended. My code also reproduces the expected result from the case without any pressure.