How can I include a finite integral in a system of non-linear 2nd order differential equation using bvp4c function and solve numerically? Tried a test with Caughey example for large-amplitude whirling of an elastic string

1 次查看（过去 30 天）
Menfin
Menfin 2020-3-28
评论： Menfin 2020-3-29
In fine, I looking for a numerical solution using bvp4c of a 3rd order differential for a clamped circular membrane where the rhs term depends on a finite integral of the soolution y.
The test example is described in the book "Solving-odes-with-matlab" from Shampine Gladwell Thompson example 1.19 et 3.9, methodology is provide but no correction is given.
Quote ------
Caughy (1970) describes the large-amplitude whirling of an elastic string by a BVP consisting of the ODE:
with boundary conditions
Here α is a physical constant with 0 <α<1. Because the whirling frequency ω is to be determined as part of solving the BVP, there must be another boundary condition. Caughey specifies the amplitude ε of the solution at the origin: µ(0) = ε.
An unusual aspect of this problem is that an important constant H is defined in terms of the solution µ(x) throughout the interval of integration:
Formulate this BVP in standard form. As in the Sturm–Liouville example, you can introduce a new variable y3(x), a first-order ODE, and a boundary condition to deal with the integral term in the definition of H. The trick to dealing with H is to let it be a new variable y4(x). It is a constant, so this new variable satisfies the first-order differential equation y4' =0. It is given the correct constant value by the boundary condition resulting from the definition of H:
Unquote ------
I tried to code several ways, but always get one od the variable y(3) or y(4) out of bound. I miss some methodlogy to get a correct output. Thanks for any help.
The following code include following analysis for differential equation:
and
y4(x)=H=cste and
and boundary conditions
and
function dydx = bvpfcn(x,y,OMEGA)
ALPHA=0.5 ; EPSILON=0.5;
dydx = zeros(4,1);
dydx = [y(2);...
-OMEGA.^2.*(((1-ALPHA.^2)./y(4)).*(1./(1+y(1).^2).^0.5)+ALPHA.^2).*y(1);...
(1./(1+y(1).^2).^0.5);...
0];
end
function res = bcfcn(ya,yb,OMEGA)
ALPHA=0.5 ; EPSILON=0.5;
res = zeros(4,1);
res = [ya(2);...
yb(2);...
yb(3)-(asinh(yb(1))-asinh(ya(1)));...
yb(4)-((1./ALPHA.^2).*(1-(1-ALPHA.^2).*yb(3)))];
end
function g = guess(x)
g = [cos(pi*x);-pi*sin(pi*x)];
end
xmesh = linspace(0,1,10);
OMEGA=pi; % initial condition
solinit = bvpinit(xmesh,@guess,OMEGA);
tol=1E-3;
options = bvpset("RelTol",tol,"AbsTol",tol,"Stats","on",'Nmax',10000);
%
tic;
Testsol4c = bvp4c(@bvpfcn, @bcfcn, solinit);
toc;

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by