Using bvp4c for coupled second order ODEs. Issue in calling bvp4c.
조회 수: 1 (최근 30일)
이전 댓글 표시
I have a system of two 2nd order ODEs. The set up of the ODE system is in the function bvp_function. I get the following error messages when I run the code:
Error using bvparguments
Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT):
The derivative function ODEFUN should return a column vector of length 100.
Error in bvp4c (line 122)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in NewPoiseuilleBaseFlow (line 38)
sol = bvp4c(@(y,X) bvp_function(y,X,b), @(Xa,Xb) boundary_conditons_fcn(Xa,Xb), guess);
Here is the code. Anyone know how to solve this problem?
% Laminar base flow for Poiseuille flow
% using Deka 2023
clear
clc
%% Setting variables
N = 100; % number of collocation points
Mach = 2.0; % Mach number
rho0 = 1.225; % density [kg/m^3]
h = 1; % height of half channel (full channel height = 2h) [m]
gamma = 1.4;
R = 287;
T0 = 288.5; % temperature (and wall temperature) (15 degC) [K]
mu0 = 1.802e-5; % dynamic viscosity [kg m^-1 s^-1]
nu0 = mu0/rho0; % kinematic viscosity [m^2/s]
kappa0 = 0.02476; % thermal conductivity [W/mK]
u0 = Mach*sqrt(gamma*R*T0); % speed [m/s]
Re = rho0*u0*2*h/mu0; % reynolds number
Pr = gamma*R*mu0/((gamma-1)*kappa0); % Prandtl number
b = Pr*(gamma-1)*Mach^2; % constant for setting up ode system
y = linspace(-1,1,N)'; % y vector
T = linspace(T0,T0,N)'; % T vector
u = linspace(1,1,N)'; % initial u vector
%% Solving the x momentum and energy equations using ode45
% initial guess for the solver
uinit = incompVelProfile(N,Mach)';
guess = bvpinit(y, uinit);
%% Solving the system of ODEs using bvp4c
sol = bvp4c(@(y,X) bvp_fn(y,X,b), @(Xa,Xb) BCs(Xa,Xb), guess);
%% Function setting up matrix system of ODEs
function dxdy= bvp_fn(y,X,b)
% X(1)=u, X(2)=u', X(3)=T, X(4)=T'
dxdy(1) = X(2);
dxdy(3) = X(4);
dxdy(4) = -0.5*X(3)*(X(4))^2-b*(X(2))^2;
dxdy(2) = -2*(X(3))^(-3/2)-0.5*(X(3))^(-1)*X(4)*X(2);
dxdy = dxdy(:);
end
%% Function for boundary conditions
function res = BCs(Xa,Xb)
% u=0, T=1 at y=-1 and u=0, T=1 at y=1.
res = [Xa(1)-0;
Xb(1)-0
Xa(3)-1
Xb(3)-1];
end
%% Function for incompressible velocity profile
function [u] = incompVelProfile(N,Mach)
% N = number of collocation points
% Mach = desired Mach number
% returns the incompressible velocity profile for channel flow
% setting variables
rho = 1.225; % density of air [kg/m-3]
nu = 1.45e-5; % kinematic viscosity [m^2/s]
T = 288; % temperature of air, 15 degC [K]
gamma = 1.4; % ratio of specific heat capacities
R = 287; % universal gas constant for air [J/(kg.K)]
% setting domain
y = linspace(-1,1,N); % y points
% working out the equivalent pressure gradient
c = sqrt(gamma*R*T); % speed of sound [m/s]
umax = Mach*c; % max velocity (centre of channel) [m\s]
dPdx = -2*umax*rho*nu; % pressure gradient required for centreline vel [Pa/m]
% obtaining velocity profile
u = zeros(N); % initialising vector
u = -0.5.*1./(rho.*nu).*dPdx.*(1-y.^2);
% have to normalise it
u = u./umax;
end
댓글 수: 0
채택된 답변
Torsten
2024년 1월 23일
guess = bvpinit(y, uinit);
y must be the spatial mesh,
uinit must either be a function handle that returns a 4x1 vector of initial values for the 4 solution components of X with an y value as input
or
a constant 4x1 vector that is assigned as initial value for the 4 solution components for each y value.
In your case, uinit seems to be an NxN matrix.
Read here about what uinit can be set to:
댓글 수: 0
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!