Hi all, I am trying to solve a fourth order differential equation using bvp4c.Below is my code.Can anyone help me to find the mistake in my code for which matlab is showing the following error?Thanks you

function sol = Fig1

solinit=bvpinit(linspace(0,30,50),@guess);

sol=bvp4c(@cantileverode,@cantileverbc,solinit);

xint=linspace(0,30,50);

Sxint=deval(sol,xint);

axis('auto')

function dy=cantileverode(x,y,l,F)

         l=17.6e-3;
         F=0.026;
         dy=zeros(4,1);
         dy(1)=y(2);
         dy(2)=y(3);
         dy(3)=y(4);
         dy(4)=(1/l^2)*(y(3)+F*cos(y(1)));
function res=cantileverbc(ya,yb,l,F)
         res=[ya(1)
             ya(2)
             yb(3)
             yb(4)];

function v =guess(x)

v=[(30-x)^4*x^2*exp(x);-4*(30-x)^3*x^2*exp(x)+(30-x)^4*(x^2*exp(x)+2*x*exp(x))];

and the errors are as follows;

Attempted to access y(3); index out of bounds because numel(y)=2.

Error in Fig1>cantileverode (line 10) dy(2)=y(3);

Error in bvparguments (line 105) testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 129) [n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] = ...

Error in Fig1 (line 3) sol=bvp4c(@cantileverode,@cantileverbc,solinit);

Answer 1

[Torsten]

You will have to supply guesses for y, y', y'' and y'''.

Your guess function only contains two components.

Best wishes

Torsten.

[Richardson Joseph]

Thank you Torsten,

I tried by providing y'' and y''',now it is giving me the following errors,

Error using vertcat Dimensions of matrices being concatenated are not consistent.

Error in Fig1>guess (line 21) v=[(30-x)^4*x^2*exp(x);

Error in bvpinit (line 104) w = feval(v,x(1),extraArgs{:});

Error in Fig1 (line 2) solinit=bvpinit(linspace(0,30,50),@guess);

[Torsten]

If you supplied a column vector of dimension (4x1) in guess, I also don't understand the error message.

Best wishes

Torsten.

[Richardson Joseph]

I tried by defining a column vector of (4x1),but still the error persist

[Torsten]

Could you again post the code you are using now ?

Best wishes

Torsten.

[Richardson Joseph]

ok,it is as below....thanks

function sol = Fig1

solinit=bvpinit(linspace(0,30,50),@guess);

sol=bvp4c(@cantileverode,@cantileverbc,solinit);

xint=linspace(0,30,50);

Sxint=deval(sol,xint);

axis('auto')

function dy=cantileverode(x,y,l,F)

        l=17.6e-3;

        F=0.026;

       dy=zeros(4,1);
         dy(1)=y(2);
         dy(2)=y(3);
         dy(3)=y(4);
         dy(4)=(1/l^2)*(y(3)+F*cos(y(1)));
function res=cantileverbc(ya,yb,l,F)
         res=[ya(1)
             ya(2)
             yb(3)
             yb(4)];
function v =guess(x)
v=zeros(4,1);
v=[(30-x)^4*x^2*exp(x);
    -4*(30-x)^3*x^2*exp(x)+(30-x)^4*(x^2*exp(x)+2*x*exp(x));
    [ 4*x^2*exp(x)*(x - 30)^3 + x^2*exp(x)*(x - 30)^4 + 2*x*exp(x)*(x - 30)^4, 24*exp(x)*(x - 30)^3 + 6*exp(x)*(x - 30)^4 + 12*x^2*exp(x)*(2*x - 60) + 36*x^2*exp(x)*(x - 30)^2 + 12*x^2*exp(x)*(x - 30)^3 + x^2*exp(x)*(x - 30)^4 + 72*x*exp(x)*(x - 30)^2 + 48*x*exp(x)*(x - 30)^3 + 6*x*exp(x)*(x - 30)^4];
    [ 4*x^2*exp(x)*(x - 30)^3 + x^2*exp(x)*(x - 30)^4 + 2*x*exp(x)*(x - 30)^4, 24*exp(x)*(x - 30)^3 + 6*exp(x)*(x - 30)^4 + 12*x^2*exp(x)*(2*x - 60) + 36*x^2*exp(x)*(x - 30)^2 + 12*x^2*exp(x)*(x - 30)^3 + x^2*exp(x)*(x - 30)^4 + 72*x*exp(x)*(x - 30)^2 + 48*x*exp(x)*(x - 30)^3 + 6*x*exp(x)*(x - 30)^4, 24*x^2*exp(x) + 144*exp(x)*(x - 30)^2 + 96*exp(x)*(x - 30)^3 + 12*exp(x)*(x - 30)^4 + 48*x^2*exp(x)*(2*x - 60) + 72*x^2*exp(x)*(x - 30)^2 + 16*x^2*exp(x)*(x - 30)^3 + x^2*exp(x)*(x - 30)^4 + 96*x*exp(x)*(2*x - 60) + 288*x*exp(x)*(x - 30)^2 + 96*x*exp(x)*(x - 30)^3 + 8*x*exp(x)*(x - 30)^4]];

Answer 2

[Torsten]

function main 
l=17.6e-3;
F=0.026;
solinit=bvpinit(linspace(0,30,50),@guess);
sol=bvp4c(@(x,y)cantileverode(x,y,l,F),@(ya,yb)cantileverbc(ya,yb,l,F),solinit);
xint=linspace(0,30,50);
Sxint=deval(sol,xint);
axis('auto')

function dy=cantileverode(x,y,l,F)
dy=[y(2)
    y(3)
    y(4)
    (1/l^2)*(y(3)+F*cos(y(1)))];

function res=cantileverbc(ya,yb,l,F)
res=[ya(1)
     ya(2)
     yb(3)
     yb(4)];

function v =guess(x)
v=[(30-x)^4*x^2*exp(x)
    -4*(30-x)^3*x^2*exp(x)+(30-x)^4*(x^2*exp(x)+2*x*exp(x))
    Insert v''
    Insert v'''];

I'm too lazy to calculate v'' and v'''.

You will have to insert them in the last two lines of the code.

Best wishes

Torsten.

[Richardson Joseph]

Thank you.It seems that now it can work,but it is showing me a warning as below;

Warning: Unable to meet the tolerance without using more than 2500 mesh points. The last mesh of 1324 points and the solution are available in the output argument. The maximum residual is 21.7711, while requested accuracy is 0.001. > In bvp4c at 323 In main at 5

Could you please give me some suggestion pertaining this warning? Thanks again

[Torsten]

Does the solution look reasonable despite the warning?

If yes: Try to increase NMax in the option settings for bvp4c:

options = bvpset('NMax',10000);
sol=bvp4c(@(x,y)cantileverode(x,y,l,F),@(ya,yb)cantileverbc(ya,yb,l,F),solinit,options);

Best wishes

Torsten.

[Richardson Joseph]

actually the figure shown is blank,I mean it only shows the axes without any plot.I am quite surprised even though I have changed the mesh points as 10000.

[Torsten]

You used a command like

xint=linspace(0,30,50);
Sxint=deval(sol,xint);
plot(xint,Sxint(1,:))

?

Best wishes

Torsten.

[Richardson Joseph]

Yes,there was a slight error in the command.Now it is showing some results. thanks alot for your kind help

[Richardson Joseph]

Hi Torsten,

I am trying to solve a very similar problem like before but now the equation is of sixth order.I wrote a code as below;

function main

l=17.6e-1;

F=0.2604;

L=5;

t=3;

solinit=bvpinit(linspace(0,5,100),@guess);

options = bvpset('NMax',10000);

sol=bvp4c(@(x,y)timoshenkoode(x,y,l,F,L,t),@(ya,yb)timoshenkobc(ya,yb,l,F,L,t),solinit,options);

xint=linspace(0,L,50);

Sxint=deval(sol,xint);

axis(['auto']);

plot(xint,Sxint(1,:))

xlabel('x')

ylabel('w')

function dy=timoshenkoode(x,y,l,F,L,t)

dy=[y(1)

    y(2)
    y(3)
    y(4)
    y(5)
    y(6)
    (1/l^2)*(F+3*t*(y(2))^2*y(3))-t*((y(3))^3+(y(2))^2*y(5)+4*y(2)*y(3)*y(4))];

function res=timoshenkobc(ya,yb,l,F,L,t)

res=[ya(1)

     ya(2)

     yb(3)

     yb(3)-l^2*yb(5)+t*0.5*l^2*(yb(2))^2*yb(3)

     yb(4)-l^2*yb(6)+t*l^2*((yb(2))^2*yb(3)+(y(2))^2*y(4))+t*0.5*(yb(2))^3

     yb(4)-t*0.5*(yb(2))^3];

 function v =guess(x)

L=5;

v=[x^2*exp(x)*(L - x)^6

    2*x*exp(x)*(L - x)^6 - 6*x^2*exp(x)*(L - x)^5 + x^2*exp(x)*(L - x)^6

    2*exp(x)*(L - x)^6 - 24*x*exp(x)*(L - x)^5 + 4*x*exp(x)*(L - x)^6 + 30*x^2*exp(x)*(L - x)^4 - 12*x^2*exp(x)*(L - x)^5 + x^2*exp(x)*(L - x)^6

    6*exp(x)*(L - x)^6 - 36*exp(x)*(L - x)^5 + 180*x*exp(x)*(L - x)^4 - 72*x*exp(x)*(L - x)^5 + 6*x*exp(x)*(L - x)^6 - 120*x^2*exp(x)*(L - x)^3 + 90*x^2*exp(x)*(L - x)^4 - 18*x^2*exp(x)*(L - x)^5 + x^2*exp(x)*(L - x)^6

    360*exp(x)*(L - x)^4 - 144*exp(x)*(L - x)^5 + 12*exp(x)*(L - x)^6 - 960*x*exp(x)*(L - x)^3 + 720*x*exp(x)*(L - x)^4 - 144*x*exp(x)*(L - x)^5 + 8*x*exp(x)*(L - x)^6 + 360*x^2*exp(x)*(L - x)^2 - 480*x^2*exp(x)*(L - x)^3 + 180*x^2*exp(x)*(L - x)^4 - 24*x^2*exp(x)*(L - x)^5 + x^2*exp(x)*(L - x)^6

    1800*exp(x)*(L - x)^4 - 2400*exp(x)*(L - x)^3 - 360*exp(x)*(L - x)^5 + 20*exp(x)*(L - x)^6 + 3600*x*exp(x)*(L - x)^2 - 4800*x*exp(x)*(L - x)^3 + 1800*x*exp(x)*(L - x)^4 - 240*x*exp(x)*(L - x)^5 + 10*x*exp(x)*(L - x)^6 - 360*x^2*exp(x)*(2*L - 2*x) + 1800*x^2*exp(x)*(L - x)^2 - 1200*x^2*exp(x)*(L - x)^3 + 300*x^2*exp(x)*(L - x)^4 - 30*x^2*exp(x)*(L - x)^5 + x^2*exp(x)*(L - x)^6];

and the error is as follows; Undefined function 'y' for input arguments of type 'double'. Error in timoshenkobc (line 2) res=[ya(1)

Error in main>@(ya,yb)timoshenkobc(ya,yb,l,F,L,t) (line 8) sol=bvp4c(@(x,y)timoshenkoode(x,y,l,F,L,t),@(ya,yb)timoshenkobc(ya,yb,l,F,L,t),solinit,options);

Error in bvparguments (line 106) testBC = bc(ya,yb,bcExtras{:});

Error in bvp4c (line 129) [n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] = ...

Error in main (line 8) sol=bvp4c(@(x,y)timoshenkoode(x,y,l,F,L,t),@(ya,yb)timoshenkobc(ya,yb,l,F,L,t),solinit,options);

[Torsten]

The two y's in the line

yb(4)-l^2*yb(6)+t*l^2*((yb(2))^2*yb(3)+(y(2))^2*y(4))+t*0.5*(yb(2))^3

must either read ya or yb.

Best wishes

Torsten.

[Richardson Joseph]

Thanks.I rectified this mistake.Now I tried to run the code and it shows me the following error;

Error using bvp4c (line 251) Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in main (line 8) sol=bvp4c(@(x,y)timoshenkoode(x,y,l,F,L,t),@(ya,yb)timoshenkobc(ya,yb,l,F,L,t),solinit,options);

[Torsten]

Your function "timoshenkoode" has 7 components whereas you only supply 6 components for boundary conditions and initial guesses.

Best wishes

Torsten.

[Richardson Joseph]

I changed the components number to "6" but still the same error of "singular jacobian".Is it related to wrong initial guess?

Thanks

[Torsten]

Usually it means that the boundary conditions are insufficient to fix a unique solution.

Best wishes

Torsten.