Hi all I am trying to plot a second order BVP in MATLAB

Question

Michael Hughes 2021-8-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/896457-hi-all-i-am-trying-to-plot-a-second-order-bvp-in-matlab

评论： darova 2021-8-12

Me and my friend are trying to plot the following BVP:

-(L)* ((d^2)y/d(x^2)) - 5*(dy/dx)/x + (a*y)/L + (c/(2*L))*(x^4)*(y^3) = 0

Assume L=6*10^-8, a=44 and c=45 (but all 3 letters can be any constant). y is a function with respect to x. The boundary conditions are y(0)=0 and y(1)=C. Where C is some constant.

Me and my friend have tried ODE45 and bvp4c but doesn't work because of (c/(2*L))*(x^4)*(y^3) and 5*(dy/dx)/x, which gives us the jacobien error.

The following below is the code we tried to use for bvp4c, we used 2 functions and one sript at the end:

Function 1:

function yprime = secondode2(x,y)
%SECONDODE: Computes the derivatives of y 1 and y 2,
%as a colum vector
a=44;  %Temperature dependent bulk constant
c=45;  %
L=6*10^-8;  %Elastic constant in the one constant approximation (splay=twist=bend)
yprime = [y(2); c/(2*L)*x^4*y(1)^3  +  a*y(1)/L  -  5*y(2)/x];
end

Function 2, for boundary conditions:

function res = boundary_conditions1(ya,yb)
res = [ya(2)
       yb(1) - 1/2];
end

Here is our sript:

guess = [1/2; 0];
xmesh = linspace(0,1,5);
solinit = bvpinit(xmesh, guess);
S = [0 0; 0 5];
options = bvpset('SingularTerm',S);
sol = bvp4c(@secondode2, @boundary_conditions1, solinit,options);
plot(sol.x,sol.y(1,:))
hold on 
title('S(r) ODE')
xlabel('x');
ylabel('solution y');