# Two couple equation first order Runge Kutta Matlab troubles?

I am solving system of two coupled equations from fluid dynamics.

$p_0 p_0'=-\dfrac{32 \beta}{R^4}$

$(p_0 p_1)'=-\dfrac{2-\sigma_v}{\sigma_v}p_0'\dfrac{8}{R}$

where $\beta$, $\sigma_v$ are constants. $R=R(z)$ radius dependent on $z$-it is not constant, where $z$ is longitudinal coordinate and $R=r-z(r-1)$, where $r$ is inlet radius. All these equations are in dimensionless form.

As a solution I need to have pressures $p_0$ and $p_1$, these pressures are first and second approximations of pressure, I got it with series expansion of pressure. And I have boundary conditions, that at the inlet $p_{0}|_{z=0}=p_i$ and at the outlet $p_0|_z=1$, for second approximation $p_1|_{z=0}=p_{1}|_{z=1}=0$. Because of this definition I implemented shooting method and got neccessary conditions.

I tried to use MATLAB code from address https://www.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html just instead of ode15 which is used on link, I need to implement ode45 (because I need to implement Runge Kutta method).

**

My problem

** is because ode45 is everywhere in matlab implemented with variables which depend on time $t$, and my pressures $p_0$ and $p_1$ are only functions on longitudinal coordinate $z$. Do I need in matlab code write my variables as $p_0(z)$ or just $p_0$, and do I need to change varible $z$ in my equations with time $t$ because matlab ode45 works only with variable $t$? Does diff in matlab accept that derivation is not according to time $t$?

I tried couple variants, and every time I got error.

Do you have any other helpful advice for solving this system of ODE with Runge Kutta method?

My code is here, I believe there is something wrong with ode45 part of the code, but I cannot conclude what.

> syms R(z) z r beta sig_v p0(z) p1(z);
>
> R(z)=r-z*(r-1);
>
> eqn1=p0(z)*diff(p0(z),z,1)==-32*beta/(R(z)^4);
> eqn2=diff(p0(z)*p1(z),z,1)==-(8*(2-sig_v)*diff(p0(z),z,1))/(sig_v*(R(z)));
> eqns=[eqn1 eqn2];
>
> vars = [p0(z); p1(z)]; origVars = length(vars);
>
> isLowIndexDAE(eqns,vars) [DAEs,DAEvars] = reduceDAEIndex(eqns,vars)
> isLowIndexDAE(DAEs,DAEvars) pDAEs = symvar(DAEs);
> pDAEvars =symvar(DAEvars);
> extraParams = setdiff(pDAEs, pDAEvars)
> %index is low, so we can go further
>
> f = daeFunction(DAEs, DAEvars, beta, r, sig_v);
> beta=1; sig_v=1; r=2;
>
> %%%%%%%%%%%% boudary condition for p and p'
> solinit = bvpinit([0 1],[0 1]);
> sol = bvp4c(@twoode,@twobc,solinit);
> xint = linspace(0,1);
> yint = deval(sol,xint);
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> % ode 45 solving
> options=odeset('RelTol', 1e-8, 'AbsTol', [1e-8 1e-8]);
> g=@(z,p0(z),p1(z), beta, R, sig_v)[eqn1; eqn2];
> [p0(z),p1(z),z]=ode45(@(z,p0(z),p1(z)) g(z,p0(z),p1(z), beta, R, sig_v), [0 1.5],[1 1/2]);

• It would help if you added the Matlab code of the latest attempt and the error message. Variable names are not important if done consistently, you can use $(t,x)$ or $(x,y)$, or $(z,p)$ or any other combination. Apr 26 '18 at 8:17

You can integrate the first equation directly, $$(p_0(z)^2)'=-\frac{64}{(2-z)^4}\implies p_0(z)^2=p_0(0)^2+\frac{64}{3(2-z)^3}-\frac83$$ As you can see, you only have one integration constant free that you can use to either fix $p_0(0)$ or $p_0(1)$, but not both at the same time.
The same goes for the integration of the second equation, you only get one integration constant that you can use to either fix $p_1(0)$ or $p_1(1)$, but it is improbable that you will get both to zero at the same time.
• My case ih that $R$ is dependent on $z$, as $p_0$ and $p_1$, so relation is $R=r-z*(r-1)$, where $r$ is constant. Also $p_0'=\dfrac{dp_0}{dz}$. How to integrate equations in this case? May 3 '18 at 13:59
• With $r=2$ you get $R=2-z$ which I inserted. And the whole computation was done with $z$ as independent variable, so nothing changes. May 3 '18 at 15:22