Solve system of two differential equations with bpv4c

I need to solve system of two equations with bvp4c method, where p is the pressure and it is only function of longitudinal coordinate z. It means p = f(z) and dpdz is the first derivative of p:

p0' = - 32 .* beta .* m0 ./ (R .^ 4 .* p0)

p1' = - ( - 8 .* p0' ./ R - p0' .* p1 - 32 .* beta .* m1 ./ R .^ 4 ) ./ p0

I defined it in Matlab function file bvpfcn.m with neccessary constants, where I want to have solution of every equation in 1001 points, it means 1x1001 is dimension of z. Exactly I need that z is going from 0 to 1 with 1001 points between 0 and 1, but I dont know where to write that:

function dpdz = bvpfcn(z,p)

beta = 1;
m0 = 1;
m1 = 0.1;
ri = 0.7;
R = ri - z .* (ri - 1);

dpdz = zeros(2,1001);
dpdz = [- 32 .* beta .* m0 ./ (R .^ 4 .* p(1))
       -( - 8 .* dpdz(1) ./ R - dpdz(1) .* p(2) - 32 .* beta .* m1 ./ R .^ 4 ) ./ p(1);
];
end

I have initial conditions where p(1) = p0 and p(2) = p1, and conditions are

p(1)|(z=0) = 8

p(1)|(z=1) = 1

p(2)|(z=0) = 0

p(2)|(z=1) = 0

which I defined if function file bcfcn.m, where pa are values of p at z=0, and pb are values of p at z=1. So boundary conditions are given as [[p0[0] p1[0]] [p0[1] p1[1]]]:

function res = bcfcn(pa,pb)
res = [[pa(1)-8    pa(1) ]
       [pb(1)-1    pb(1) ]];
end

My solution need to be of shape:

Because of that I am making guess function of cosine type, I am not sure is this good assumption?

function g = guess(z)
g = [cos(z)
     cos(z)];
end

And call all with code:

clear all;
beta = 1;
m0 = 1;
m1 = 0.1;
ri = 0.7;
xmesh = linspace(0,pi/2,1001);
solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@bvpfcn, @bcfcn, solinit);

This is error which I got:

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in bvp4c>colloc_RHS (line 600)
Phi(1:nBCs) = Gbc(y(:,Lidx),y(:,Ridx),ExtraArgs{1:nExtraArgs});

Error in bvp4c (line 189)
  [RHS,yp,Fmid,NF] = colloc_RHS(n,x,Y,ode,bc,npar,xyVectorized,mbcidx,nExtraArgs,ExtraArgs);

Answer 1

[darova]

First of all: it's a bad practice to assign variables

dpdz = [- 32 .* beta .* m0 ./ (R .^ 4 .* p(1))
       -( - 8 .* dpdz(1) ./ R - dpdz(1) .* p(2) - 32 .* beta .* m1 ./ R .^ 4 ) ./ p(1);
];

Example

>> [1 -2]
ans =
     1    -2
>> [1 - 2]
ans =
    -1

Try (no need of element-wise operators .* and ./)

dpdz(1,1) = -32*beta*m0/R^4*p(1);
dpdz(2,1) = (8*dpdz(1)/R + dpdz(1)*p(2) + 32*beta*m1/R^4)/p(1);

Boundary conditions should be a column

function res = bcfcn(pa,pb)
res = pa(1)-8    
        pa(1)
        pb(1)-1    
        pb(1)];
end

Your equations look like usual ode. Why are they bvp?

I tried this

xmesh = linspace(0,0.01,1001);
[z,p] = ode45(@bvpfcn,xmesh,[3 3]);
plot(z,p)

result

[I G]

They are bvp, because I have two boundary conditions, it means for every value of p I have assigned vaules (boundary conditions) at coordinate z=0 and at z=1. Also solution which you get with ode is not ok, because I need to get plot as given in my first question. Also now with your version of boundary conditions I got error:

Error using bvparguments (line 111)
Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT):
  The boundary condition function BCFUN should return a column vector of length 2.

Error in bvp4c (line 130)
    bvparguments(solver_name,ode,bc,solinit,options,varargin);

I dont know how to make column vector of length 2 with 4 bounday conditions?

[darova]

• I dont know how to make column vector of length 2 with 4 bounday conditions?

I don't know too!

If you have first order ODE

p0' = - 32 .* beta .* m0 ./ (R .^ 4 .* p0)

You can't change end boundary if you have one initial

p0(i+1) = p0(i) + p0'*dt

Something is missing