Prev: Line plot starting at second x axis tick and ...
Next: How to call a program after figure becomes visible
From: Dave on 5 Jun 2010 09:53
Dear all,

I am trying to solve a system of ODE's with boundary conditions using bvp5c. I was wondering if anybody knows how to impose two sets of boundary conditions. I need

y(0)=13.4171, y(2*pi)=13.4171 for y(1)

and

y(0)=115, y(2*pi)=115 for y(2)

I get an error stating "The boundary condition function BCFUN should return a column vector of length 2."

function BVPattempt

r=10; % radius of circular flow (kpc)
alpha=0.11667; % angle of spiral inclination (radians)

solinit = bvpinit(linspace(0,alpha*r*pi,1000),[13.4717 115]);

sol = bvp5c(@f,@ex1bc,solinit);


% The solution at the mesh points
x = sol.x;
y = sol.y;


figure;
plot(x,y(2,:)')
hold on
%plot(x,y(2,:)')
legend('u','v');
xlabel('eta');
ylabel('u/v');
title('None');
%---------------------------------------------------------------

function dudy = f(y,u)
dudy = zeros(2,1);

% define physical constants

omega=25; % angular velocity at omega(r) (km /sec /kpc)
kappa=31.3; % epicyclic frequency (km /sec /kpc)
omegapattern=13.5; % angular pattern velocity (km /sec /kpc)
c=8.56; % speed of sound (km /sec)
r=10; % radius of circular flow (kpc)
alpha=0.11667; % angle of spiral inclination (radians)
A=72.92; % amplitude (km /sec)^2

v0=r*(omega - omegapattern);
u0=alpha*r*(omega - omegapattern);

dudy(1) = (u(1)/((u(1)^2) - (c^2)))*(2*omega*(u(2) - v0) + ((2*A)/(alpha*r))*sin(((2*y)/(alpha*r))));
dudy(2) = -((kappa^2)/(2*omega))*((u(1) - u0)/u(1));


%---------------------------------------------------------------

% Not sure how to impose two sets of boundary conditions

function res = ex1bc(ya,yb)

res = [ ya(1) - 13.4171
yb(1) - 13.4171
ya(2) - 115
yb(2) - 115];

Also there is a discontinuity which the scheme seems to struggle with. Can anybody recommend how to deal with discontinuities please? I have looked at the shockbvp demo code but it isn't much help
From: Torsten Hennig on 6 Jun 2010 22:37
> Dear all,
>
> I am trying to solve a system of ODE's with boundary
> conditions using bvp5c. I was wondering if anybody
> knows how to impose two sets of boundary conditions.
> I need
>
> y(0)=13.4171, y(2*pi)=13.4171 for y(1)
>
> and
>
> y(0)=115, y(2*pi)=115 for y(2)
>
> I get an error stating "The boundary condition
> function BCFUN should return a column vector of
> length 2."
>
> function BVPattempt
>
> r=10; % radius of circular flow (kpc)
> alpha=0.11667; % angle of spiral inclination
> (radians)
>
> solinit =
> bvpinit(linspace(0,alpha*r*pi,1000),[13.4717 115]);
>
> sol = bvp5c(@f,@ex1bc,solinit);
>
>
> % The solution at the mesh points
> x = sol.x;
> y = sol.y;
>
>
> figure;
> plot(x,y(2,:)')
> hold on
> %plot(x,y(2,:)')
> legend('u','v');
> xlabel('eta');
> ylabel('u/v');
> title('None');
> %-----------------------------------------------------
> ----------
>
> function dudy = f(y,u)
> dudy = zeros(2,1);
>
> % define physical constants
>
> omega=25; % angular velocity at omega(r)
> (km /sec /kpc)
> kappa=31.3; % epicyclic frequency (km /sec
> /kpc)
> omegapattern=13.5; % angular pattern velocity (km
> /sec /kpc)
> c=8.56; % speed of sound (km /sec)
> r=10; % radius of circular flow (kpc)
> alpha=0.11667; % angle of spiral inclination
> (radians)
> A=72.92; % amplitude (km /sec)^2
>
> v0=r*(omega - omegapattern);
> u0=alpha*r*(omega - omegapattern);
>
> dudy(1) = (u(1)/((u(1)^2) - (c^2)))*(2*omega*(u(2) -
> v0) + ((2*A)/(alpha*r))*sin(((2*y)/(alpha*r))));
> dudy(2) = -((kappa^2)/(2*omega))*((u(1) - u0)/u(1));
>
>
> %-----------------------------------------------------
> ----------
>
> % Not sure how to impose two sets of boundary
> conditions
>
> function res = ex1bc(ya,yb)
>
> res = [ ya(1) - 13.4171
> yb(1) - 13.4171
> ya(2) - 115
> yb(2) - 115];
>
> Also there is a discontinuity which the scheme seems
> to struggle with. Can anybody recommend how to deal
> with discontinuities please? I have looked at the
> shockbvp demo code but it isn't much help

Two first-order ODEs need _two_ boundary conditions,
not four.
If the solution comes out to be periodic with period
2*pi, the boundary conditions at 2*pi will automatically
be satisfied.

Best wishes
Torsten.
 | 
Pages: 1
Prev: Line plot starting at second x axis tick and ...
Next: How to call a program after figure becomes visible