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]);
  • 2
    $\begingroup$ 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. $\endgroup$ 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.

Usually a boundary value problem requires an order 2 ODE, then you would get 2 integration constants in every step that you can (try to) solve for the 2 boundary conditions.

  • $\begingroup$ 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? $\endgroup$
    – nick_name
    May 3 '18 at 13:59
  • 1
    $\begingroup$ 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. $\endgroup$ May 3 '18 at 15:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.