solving 4th order ode

Hello everyone! I would appreciate your help in this section.

I have this 4th order,non-linear ,homogeneous ode that I would like to solve:

(2B+1+y)y''-k(1+y)y''''=0 with: y(0)=y(L)=0;y'(0)=y'(L)=1; where B,k and L are constants.

This is my matlab code so far:

L=1;B=0.25;k=0.1;
dydt = zeros(4,1);
dydt(1)=y(2);
dydt(2)=y(3);
dydt(3)=y(4);
dydt(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));
y1(1)=0; y1(end)=0; y2(1)=1; y2(end)=1;yint=[y1(1);y1(end);y2(1);y2(end)];
tspan = [0 1];
[t,y] = ode45(@(t,y) dydt, tspan, yint);
plot(t,y,'-o')

I think that something is not right. I have almost no expericance with solving ode's with MATLAB, especially with this high order so basically I gathered pieces of codes from things I saw online with the hope that it will assemble to a reasonable solution. I'm not sure I'm in a good path at all. I would be grateful for some guidance. Roi

Answer 1

[darova]

Use bvp4c()

You already have written your equation as system of first-order equations

function dydt = ode4(t,y)
    L=1;B=0.25;k=0.1;
    dydt = zeros(4,1);
    dydt(1)=y(2);
    dydt(2)=y(3);
    dydt(3)=y(4);
    dydt(4)=(2.*B+1+y(1)).*y(3))./(k.*(1+y(1));
end

Here is function for your boundary conditions

function res = bc4(y0, y1)
    res = [y0(1)-0        % y(0) = 0
           y0(2)-1        % y'(0) = 1
           y1(1)-0        % y(end) = 0
           y1(2)-1];      % y'(end) = 1
end

See bvp4c

[Roi Bar-On]

Hey Darova and thank you so much for your quick answer!

This is the updated code:

function dydt = bvpfcn(x,y)

dydt = zeros(4,1);

dydt(1)=y(2);

dydt(2)=y(3);

dydt(3)=y(4);

%L=1;

B=0.25;k=0.1;

dydt(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));

end

function res = bcfcn(y0, y1)

res = [y0(1)-0 % y(0) = 0

y0(2)-1 % y'(0) = 1

y1(1)-0 % y(end) = 0

y1(2)-1]; % y'(end) = 1

end

Now I'm getting :

Not enough input arguments.

Error in bvpfcn (line 5)

dydt(1)=y(2);

I think that because it's a 4th order ode then maybe bvpfcn(x,y) shold be a function of y's derivatives as well? I'm really not sure.

[darova]

Where is the main code?

0_0

[Roi Bar-On]

Here it is:

function dydx = bvpfcn01(x,y) % equation to solve

dydx = zeros(4,1);

dydx(1)=y(2);

dydx(2)=y(3);

dydx(3)=y(4);

%L=1;

B=0.25;k=0.1;

dydx(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));

end

%--------------------------------

function res = bcfcn(y0, y1)

res = [y0(1)-0 % y(0) = 0

y0(2)-1 % y'(0) = 1

y1(1)-0 % y(end) = 0

y1(2)-1]; % y'(end) = 1

end

%--------------------------------

function g = guess(x) % initial guess for y and y'

g = [x.^4

x.^3

x.^2

x];

xmesh = linspace(0,0.01,1);

solinit = bvpinit(xmesh, @guess);

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

plot(sol.x, sol.y, '-o')

end

%--------------------------------

%--------------------------------

I still don't know that to take as the arguments of bvpfcn(x,y). are those guesses of x and y? As for the g function, I just took arbitrary x^n.

[darova]

Try this:

function main
    xmesh = linspace(0,1,50);
    solinit = bvpinit(xmesh, [0 0 0 0]);
    sol = bvp4c(@bvpfcn01, @bcfcn, solinit);
    % y1 = sol.y(3,:);
    plot(sol.x, sol.y, '-o')
end

And please use button for code inserting:

[Roi Bar-On]

I think I'm missing something. This is the hall code:

function dydx = bvpfcn01(y) % equation to solve

dydx = zeros(4,1);

dydx(1)=y(2);

dydx(2)=y(3);

dydx(3)=y(4);

%L=1;

B=0.25;k=0.1;

dydx(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));

end

%--------------------------------

function res = bcfcn(y0, y1)

res = [y0(1)-0 % y(0) = 0

y0(2)-1 % y'(0) = 1

y1(1)-0 % y(end) = 0

y1(2)-1]; % y'(end) = 1

end

%--------------------------------

function main

xmesh = linspace(0,1,50);

solinit = bvpinit(xmesh, [0 0 0 0]);

sol = bvp4c(@bvpfcn01, @bcfcn, solinit);

% y1 = sol.y(3,:);

plot(sol.x, sol.y, '-o')

end

%--------------------------------

I'm getting this error:

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

Error in bvpfcn01 (line 4)

dydx(1)=y(2);

From what I understand MATLAB can't find y(2) for instance because ''y'' is not defined as an array, therefore the numel comment. How can I define it? and when I type in the command window bvpfcn01(y) , should I submit a value of y as an initial guess? Because it doesn't work. Shouldn't it be bvpfcn01(xmesh,y)? I'm really clueless here:( .. Thanks for your patience.

[darova]

Should be two arguments, even if x is not used

function dydx = bvpfcn01(x,y) % equation to solve
% OR
function dydx = bvpfcn01(~,y) % equation to solve

[Roi Bar-On]

Okay but do I need to supply an initial guess for both of them? Because they're not defined.... Also this numel issue with the array of y.

[Roi Bar-On]

I changed the code a bit and now it works!

This is the new code after arragment:

function main
    xmesh = linspace(0,1,50);
    solinit = bvpinit(xmesh, [0 0 0 0]);
    sol = bvp4c(@bvpfcn01, @bcfcn, solinit);
    % y1 = sol.y(3,:);
    plot(sol.x, sol.y(3,:), '-o')
end
%--------------------------------


function dydx = bvpfcn01(~,y) % equation to solve
dydx = zeros(4,1);
    dydx(1)=y(2);
    dydx(2)=y(3);
    dydx(3)=y(4);
    %L=1;
    B=0.25;k=0.1;
    dydx(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));
end
%--------------------------------
function res = bcfcn(y0, y1)
    res = [y0(1)-0        % y(0) = 0
           y0(2)-1        % y'(0) = 1
           y1(1)-0        % y(end) = 0
           y1(2)-1];      % y'(end) = 1
end

How did you know that y1 = sol.y(3,:) is the dependant variable that I'm looking for? (the original y). I figured that MATLAB gives the solution as a structure. How can I get as an output as well this structure and extract the output array of ''y'' from it?

Again, thank you so much for your help!

[darova]

Also you can save your functions (bvpfcn01 and bcfcn) as separate files

You are already extracting values using command:

% y1 = sol.y(3,:);

What do you mean?

[Roi Bar-On]

Sorry for the delay.

I meant that sol.y is the array that holds the derivatives of y. basically it contains y,y',y'',y''',y''''.

sol.y(1,:) shouldn't be y? why sol.y(3,:) is y according to what you say?

Roi

[darova]

• sol.y(1,:) shouldn't be y?

Exactly!

• sol.y(1,:) - y
• sol.y(2,:) - y'
• sol.y(3,:) - y''
• sol.y(4,:) - y'''

[Roi Bar-On]

Hi again!

Iv'e noticed that the solution that I'm getting can be only constant - because the initial guess is a scalar and not a variable. I wrote a new code but I was forced to write it down as a bvp with y(0)=0 and y(1)=1, and actually I don't know what's the value of y(1) so I marked that as ''x''. Actually I only know (as presented in the previous code) that y'(0)=y'(1)=0. This is my code:

function shooting_method
clc
clear all
x=0.5;
x1=fzero(@solver,x);
end

function F=solver(x)
options=odeset('RelTol',1e-8,'AbsTol',[1e-8,1e-8]);
[t,u]=ode45(@equation,[0 1],[0 x],options);%ode15s,ode23s are options as well
s=length(t);
F=u(s,1)-1;
figure(1)
plot(t,u(:,1))
xlabel('x')
    ylabel('rho')
    title('S.M')
end

function dy=equation(~,y)
     dy=zeros(2,1);
     dy(1)=y(2);
     A=1;E=1;
     dy(2)=1./A.*log10(E.*y(1));
    
end

I've simplified the problem now and got a more basical ode as can be seen. I marked y(1)=x , because I don't this value. I'm getting weird results. I wanted to know how in this code I can incorporate my conditions y'(0)=y'(1)=0 and not the bvp ones. I just don't know how to use bcfcn here.

I would love to get some help.

Thank you,

Roi

[darova]

Change your initial condition (maybe 0.1 or 0.01)

dy(2)=1./A.*log10(E.*y(1))

because

>> log10(0)
ans =
  -Inf

And try

x1=fsolve(@solver,x);

[Roi Bar-On]

Hi.

Unfortunately it doesn't work. It mentions that I need a certain ''toolbox''.

I didn't realize yet how do I blend in y'(0)=y'(1)=0.

Roi

[darova]

Change initial condition

x=1;
x1=fzero(@solver,x);

[Roi Bar-On]

yeah but that still holds for y(0)=0, y(1)=x. Those are Dirichlet b.c . I have Neumann b.c y'(0)=y'(1)=0 . I don't understand how to put them in the code instead of y(0)=0, y(1)=x.

[darova]

THe same way you wrote it in solver

y'(0) = 0; y(0) = ? (initial conditions)

[t,u]=ode45(@equation,[0 1],[x 0],options);%ode15s,ode23s are options as well
s=length(t);
F=u(s,2)-0;        % y'(1) = 0

[Roi Bar-On]

I didn't understand.. here in [t,u]=ode45(@equation,[0 1],[x 0],options) you just said that y(1)=0 . and what F=u(s,2)-0 means? 2 is for the derivative of y? Basically as I said I don't know that's y(0) nor y(1).. I only know y'(0)=y'(1)=0.

Thank you for your patience,

Roi

[darova]

Why don't you understand me? Using ode45 you can set only initial conditions

The using fsolve you are adding boundary condition (dy(1)=0 in your case)

Or you can use bvp4c method to solve this problem

[Roi Bar-On]

Finally I understood your explanation and the rational, sorry for being so ignorent in this topic.

I've manually ran the code step by step and I see that basically nothing really happens. Basically we're giving the initial guess x=1 and the plot is just getting values of 1 eveywhere because 1 is really a solution. It seems that the code doesn't look for any other solutions and therefore I'm not getting any physical results. This is the final code for now:

function shooting_method
clc
clear all
x=1;
x1=fzero(@solver,x);
end

function F=solver(x)
options=odeset('RelTol',1e-8,'AbsTol',[1e-8,1e-8]);
[t,u]=ode45(@equation,[0 1],[x 1e-8],options);%1e-8 is for the function and x is for the derivative%ode15s,ode23s are options as well
s=length(t);
F=u(s,2)-0; %y'(1)=0
figure(1)
plot(t,u(:,1))
xlabel('x')
    ylabel('rho')
    title('S.M')
end

function dy=equation(~,y)
     dy=zeros(2,1);
     dy(1)=y(2);
     A=1;E=1;
     dy(2)=1./A.*log10(E.*y(1));
    
end

I would be glad for some help. Maybe this method isn't stable enough or the problem is stiff? If so do you happen to know any other method that might help me solve that problem?

Thanks,

Roi

[darova]

I tried this

x=0.5;
x1=fsolve(@solver,x);

And set limits for y axis:

ylim([0 2])

THe result (y'(0)=y'(1)=0.):

Answer 2

[Roi Bar-On]

Thanks,

Roi

[Elayarani M]

Hi. I want to solve the following equation using ode45 solver.

How to set the initial conditions?