2nd order non linear ODE

I cannot figure out how to solve the folowing non linear eigen value ODE numerically on MATLAB.

I know the value of the second term, g as well. My boundary conditions are and . I know how to solve it by hand but MATLAB is a different story for me. Please help out.

[KSSV]

Have a look on ode45

[Ahsan Mujtaba]

Do I have to convert it into system of 2 first order equation?

Answer 1

[Star Strider]

‘Do I have to convert it into system of 2 first order equation?’

Yes.

The Symbolic Math Toolbox can make that easier. Use the odeToVectorField and matlabFunction functions to create it as an anonymous function.

Also, since this is a boundary value problem, bvp4c or similar functions may be most appropriate. (I have no idea what ‘E’ and ω are, however using this approach and using random values for those and guessing 9.81 for ‘g’, I got it to work with bvp4c.)

I would post my solution, however I would not want to deprive you of the experience of discovering how to do all this for yourself.

[Ahsan Mujtaba]

syms y(x)
a=0.4871;%value of E/w
g=10;

eqn=diff(y,2)==(a-(sin(sqrt(g)*x)/g))*y;%converted equation to x and y for ease
v = odeToVectorField(eqn)
M = matlabFunction(v,'vars',{'x','Y'})
interval=[-pi/2 pi/2];%sqrt(g)x interval
yInit=[0 0];
ySol = ode45(M,interval,yInit);
xValues=linspace(-pi/2,pi/2,10000);
yValues = deval(ySol,xValues,1);
plot(xValues,yValues)

The problem is I have to plot I don't think my ouptut is coming right because it is giving me the straight line the output should be sin type wave. I think I am not applying the limits rightly please help out.

[Star Strider]

With zero initial conditions for both differential equations, the integrated result will be identically zero.

Ityy is possible to get a non-zero result with non-zero initial conditions:

yInit=[0 0]+eps;

That aside, something like this could do what you want:

syms y(x)
a=0.4871;%value of E/w
g=10;
eqn=diff(y,2)==(a-(sin(sqrt(g)*x)/g))*y;%converted equation to x and y for ease
v = odeToVectorField(eqn)
M = matlabFunction(v,'vars',{'x','Y'})
interval=[-pi/2 pi/2];%sqrt(g)x interval
interval = linspace(-pi/2,pi/2,10000)*sqrt(g);
yInit=[0 0]+eps;
[x,y] = ode45(M,interval,yInit);
xValues=linspace(-pi/2,pi/2,10000)*sqrt(g);
yValues = y(:,1);
plot(xValues,yValues)

Experiment to get the result that you want.

[Ahsan Mujtaba]

with non zero initial conditions I get the results but again they are not correct. I have two boundary conditions am I solving it with the wrong function?

[Star Strider]

I forgot about it being a boundary value problem. I solved it as such a few days ago. Use bvp4c to integrate it. I used the symbolic code similar to what you are using to create the anonymous function for the differential equation system, then hand-coded the rest.

If you have problems, I can use my solution to guide you through solving it. The plot looks like a parabola (or part of a sine curve).

[Ahsan Mujtaba]

yes the plot should look like a part of sine curve please guide as well.

Many thanks

[Star Strider]

My pleasure!

See the bvp4c documentation that I linked to. It has everything you need to solve it (and is what I used, since I do not code boundary value problems very often, and neded to refer to it).

[Ahsan Mujtaba]

are you getting two plots? And are your plots changing with the value of g e.g. for g=10 there is a large parabola with endpoints at pi/2s.

[Ahsan Mujtaba]

function dydx = mat4ode(x,y,lambda) % equation being solved
g=0.1;
q = 4/(g^2);
dydx = [y(2)
      -(2*lambda - ((sin(sqrt(g)*x))^2/g))*y(1)];
end
function res = mat4bc(ya,yb,lambda) % boundary conditions
res = [ya(-4.96)
       yb(4.96)
       ya(2)-1];
end
function yinit = mat4init(x) % initial guess function it is confussing for me
yinit = [cos(4*x)
        -4*sin(4*x)];
end
lambda=1.4390; %E/w value 

solinit = bvpinit(linspace(-pi/2,pi/2,1000),@mat4init,lambda);
sol = bvp4c(@mat4ode, @mat4bc, solinit);
xint = linspace(-pi/2,pi/2,1000);
Sxint = deval(sol,xint);
plot(xint,Sxint)

title('Eigenfunction of Mathieu''s Equation.') 
xlabel('x')
ylabel('y')
legend('y','y''')

There are various things consfusing for me and this code only runs for ya(2),yb(2),ya(1)-1 BC's. How can I make it run for my BC's? Also guesses bother me. Please look at the code as well and guide.

[Star Strider]

I selected only the first row of the solution to plot, so only , with this result —

The second row simply plots a straight line with a negative slope, so I am not plotting it.

I am using my own XTickLabel strings for clarity.

[Ahsan Mujtaba]

Can you please see my code and guide me the guesses?

[Star Strider]

You must have posted your code while I was posting my plot. I did not see it before.

There are only two boundary conditions. You are coding three of them. Also, when I try to run your code, the first boundary condition throws this error:

Array indices must be positive integers or logical values.

Error in ... mat4bc (line ##)
res = [ya(-4.96)

Please understand that these are actually subscripts into the differential equation function variables, not arguments. You need to re-code the boundary conditions.

I coded only 25 points in the initial mesh. Fewer than 1000 of them may be desirable, unless your code throws an error about a singular Jacobian (that I doubt it will), in which situation more points will be necessary. The bvp4c function will tell you if it needs a larger mesh.

[Ahsan Mujtaba]

how can I code my boundary conditions I am sorry I am not very good in it.

[Star Strider]

The documentation for bvp4c describes exactly how to code them.

[Ahsan Mujtaba]

function dydx = bvpfcn(x,y)

dydx = zeros(2,1);
dydx = [y(2)
       (2*0.4878-(((sin(sqrt(10))*x)^2)/10))*y(1)];
end
function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)];
end
function g = guess(x)
g = [sin(x)
     4*cos(x)];
end
xmesh = linspace(-pi/2,pi/2,100);

solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@bvpfcn, @bcfcn, solinit);
plot(sol.x, sol.y)

Can you please see this code this time I have carefuly coded the BCs but still my plots are not coming right. I am a beginner I would appreciate any leads. Thanks

[Star Strider]

I am not certain what it is supposed to look like, or if my plot is correct. (It would appear to be correct, however I am not certain it is.)

I used a different ‘guess’, and when I use the same in your code, I get essentially the same result I got in mine, although your ‘dydx’ is different from the function I used, that being the result of the matlabFunction call in my code. I also used only 25 points in ‘xmesh’, and that made for a smoother plotted result.

UPDATE — I did not previously realise it, however Mathieu’s equation is actually in the MATLAB documentation. (I just now discovered that in an Interweb search.) See Solve BVP with Unknown Parameter for details. It woulld appear that my result is correct, given that the parameters may not be the same (I have not checked them, since it is late here, and the end of my day.).

I find that comforting!