0

I need to solve the following

-cos(y)y''+sin(y)y'^2+sin(y)=0, y'(0)=y'(1)=0, such that y=y(t)

I find it hard to solve because of the term y'^2 and also the boundary conditions.

  • What kind of solution do you hope to find? I'm not aware of a BVP solver in scilab, but you could implement a shooting method. For a symbolic solution compare with `u=f(y)` where then `u'=f'(y)y'` and `u''=f'(y)y''+f''(y)y'^2` and try to identify a suitable `f`. – Lutz Lehmann Mar 09 '21 at 12:46
  • 1
    Scilab has `bvode` and `bvodeS` to solve boundary value problems. – Stéphane Mottelet Mar 09 '21 at 13:17
  • @LutzLehmann I can do `z=(-cos(y)y')` and then `z'=-sin(y)` ? What about the boundary condition? – Hariz Khaled Mar 09 '21 at 13:30
  • With the sign change you can now simply use `u=sin(y)`, then `-u''+u=0`, with only the trivial solution. – Lutz Lehmann Mar 09 '21 at 13:36
  • @LutzLehmann Thanks! What about the general for instance, by replacing `cos(y)` by `(1-cos(y))`, I think a shooting method taking the boundary `y(0)` and `y(1)` but not for `y'(0)` and `y'(1)`? – Hariz Khaled Mar 09 '21 at 14:06
  • That would be a different question, and in its current form quite off-topic here. For the mathematical possibility you should ask in math.SE, for algorithmic problems in scicomp.SE, and only if you have problems with concrete code (after exhausting the obvious debugging steps) you should ask here. // Single-shooting works for any correctly formulated BVP, but it might be a shot-in-the-dark, too far away from the target, possibly running into singularities that the solution does not have. For that one has to make the system more Lipschitz or use multiple shooting. – Lutz Lehmann Mar 09 '21 at 14:52

1 Answers1

3

Here is the Scilab code for your bvp

-cos(y)y''+sin(y)y'^2+sin(y)=0, y'(0)=y'(1)=0, y(0)=0, y(1)=1.5

but with different boundary condition not giving the trivial solution. You have first to write y'' as a function of y and y'. The function fsub computes y'' as a function of u=[y,y']

function ysecond=fsub(x,u)
    y=u(1);
    yprime=u(2);
    ysecond = sin(y)/cos(y)*(1+yprime^2);
end

function g=gsub(i, u),
    y=u(1);    
    select i
      case 1 then  // x=zeta(1)=0
        g = y // y(0)=0
      case 2 then // x=zeta(2)=1
        g = y-1.5 // y(1)=1.5
    end
end

N=1;
m=2;
x_low=0
x_up=1;
xpoints=linspace(0,1,100);
zeta=[0,1];


u = bvodeS(xpoints,m,N,x_low,x_up,fsub,gsub,zeta)
plot(xpoints,u(1,:))

enter image description here

Stéphane Mottelet
  • 2,853
  • 1
  • 9
  • 26