Asked by migurel figueroa
on 13 Mar 2019

I can't understand bvp4c. I have a system of 432 differential equations, and obviously I have 432 dependent variables y. I know my boundary values are y(t=0) = 0 and y (t=tf) = y_f. But bvp4c asks me for that handle that includes the boundary conditions and I don't know how to put it I know it should be a bc((a)(b)) = 0. Following the example from mathworks DOESN'T HELP. my final boundary condition should be an array of 432 different values. I'm getting this error: "The boundary condition function BCFUN should return a column vector of length 432."

Here's my code

function Dynamic_condiciones_de_frontera

clear all

clc

F=Op_data_dinamico.F;

z=Op_data_dinamico.z;

V=Op_data_dinamico.V;

L=Op_data_dinamico.L;

K= Op_data_dinamico.K;

W=Op_data_dinamico.W;

U=Op_data_dinamico.U;

x=Op_data_dinamico.x;

M=Op_data_dinamico.M;

N=Op_data_dinamico.N;

NC=Op_data_dinamico.NC;

rr=2;

solinit = bvpinit (linspace(0,3600,10),[ones(1,432)]);

x_final=Op_data_dinamico.x(:);

size(x_f)

sol = bvp4c(@ex1ode, @bcval, solinit);

end

function dxdt = ex1ode(t,x)

%% Platos intermedios %%

for j=2:N-1

for i=1:NC

dxdt(j,i)=(((F(j)*z(i))+(L(j-1)*x(j-1,i))+(V(j+1)*K(j+1,i)*x(j+1,i))-((L(j)+U(j))*x(j,i))-((V(j)+W(j))*K(j,i)*x(j,i))))/M;

end

end

%% Tope %%

L_destill=L(1)/rr;

for i=1:NC

dxdt(1,i)=((V(2)*K(2,i)*x(2,i)-L(1)*x(1,i)-L_destill*x(1,i)))/M;

end

%% Fondo %%

for i=1:NC

dxdt(N,i)=(L(N-1)*x(N-1,i)-V(N)*K(N,i)*x(N,i)-L(N)*x(N,i))/M;

end

dxdt=dxdt(:);

end

function residual = bcval(xa,xb)

x_f=Op_data_dinamico.x(:);

residual = [xa(1); xb(1)-x_final];

end

Answer by Torsten
on 13 Mar 2019

As initial condition, you give a row vector of ones of length 432.

So in ex1ode, bvp4c will give you a vector of length 432 for x. So why do you expect you can work with x in ex1ode as if it were a matrix of size (N,NC) ?

And in bcval, x_final (which does not even exist because you named it x_f) must be a column vector of length 431 for that at least technically, bvp4c will be satisfied. But note that your boundary condition does not make sense at all because you want the first ODE take a value of 0 at t=0 and simultaneously all the values of x_final at t=3600.

migurel figueroa
on 13 Mar 2019

Here's how the code looks now:

What errors does it still have?

function Dynamic_condiciones_de_frontera

clear all

clc

F=Op_data_dinamico.F;

z=Op_data_dinamico.z;

V=Op_data_dinamico.V;

L=Op_data_dinamico.L;

K= Op_data_dinamico.K;

W=Op_data_dinamico.W;

U=Op_data_dinamico.U;

x=Op_data_dinamico.x;

M=Op_data_dinamico.M;

N=Op_data_dinamico.N;

NC=Op_data_dinamico.NC;

rr=2;

zinit = bvpinit (linspace(0,3600,50),[zeros(1,216) ones(1,216)]);

x_f=Op_data_dinamico.x(:);

size(x_f)

sol = bvp4c(@Dynami_comp_robusto, @bcval, zinit);

end

function dxdt = ex1ode(t,x)

%% Platos intermedios %%

for j=2:N-1

for i=1:NC

dxdt(j,i)=(((F(j)*z(i))+(L(j-1)*x(j-1,i))+(V(j+1)*K(j+1,i)*x(j+1,i))-((L(j)+U(j))*x(j,i))-((V(j)+W(j))*K(j,i)*x(j,i))))/M;

end

end

%% Tope %%

L_destill=L(1)/rr;

for i=1:NC

dxdt(1,i)=((V(2)*K(2,i)*x(2,i)-L(1)*x(1,i)-L_destill*x(1,i)))/M;

end

%% Fondo %%

for i=1:NC

dxdt(N,i)=(L(N-1)*x(N-1,i)-V(N)*K(N,i)*x(N,i)-L(N)*x(N,i))/M;

end

dxdt=dxdt(:);

end

function res = bcval(xa,xb)

x_f=Op_data_dinamico.x(:);

res= [xa(1:216);

xb(1:216)-x_f(1:216)];

end

Torsten
on 13 Mar 2019

- Why do you overwrite the xa and xb values passed in "bcval" by zero values ? That does not make sense.
- x - as passed to ex1ode - is a vector of length 432. You use it as if it were a matrix of size (N,NC). Note that all the variables and arrays taken from "Op_data_dinamico" in function "Dynamic_condiciones_de_frontera" are not visible in function "ex1ode".

migurel figueroa
on 13 Mar 2019

1- Sorry, that was an act of desperation and I already deleted it.

2- How should I write x to act as the boundary conditions for my xb from my equation 216 to 432 (I assume it must be from 216 from 432 because the initial boundary uses my first 216 equations). I managed it to give me a result, an awful result but still a result. This is how I finally ended with the boundary function:

function res = bcval(xa,xb)

x_f=Op_data_dinamico.x(:);

res= [xa(1:216);

xb(217:432)-x_f(217:432)];

end

However, it's given me this graph that doesn't make sense at all:

The values of X across time should always be between [0 and 1] and as you can see it's giving me negative numbers. What am I doing wrong?

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.