Solve system of two differential equations with bpv4c

1 visualizzazione (ultimi 30 giorni)
I G
I G il 20 Mar 2020
Commentato: darova il 23 Mar 2020
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);

Risposte (1)

darova
darova il 20 Mar 2020
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
  2 Commenti
I G
I G il 23 Mar 2020
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
darova il 23 Mar 2020
  • 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

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by