You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
solving 4th order ode
52 views (last 30 days)
Show older comments
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
Accepted Answer
darova
on 1 Oct 2019
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
21 Comments
Roi Bar-On
on 1 Oct 2019
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.
Roi Bar-On
on 1 Oct 2019
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.
Roi Bar-On
on 2 Oct 2019
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.
Roi Bar-On
on 2 Oct 2019
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
on 2 Oct 2019
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!
Roi Bar-On
on 7 Nov 2019
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
Roi Bar-On
on 8 Nov 2019
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
Roi Bar-On
on 9 Nov 2019
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
Roi Bar-On
on 10 Nov 2019
Roi Bar-On
on 10 Nov 2019
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
Roi Bar-On
on 19 Nov 2019
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
More Answers (1)
Roi Bar-On
on 7 Nov 2019
Thanks,
Roi
1 Comment
Elayarani M
on 23 Sep 2020
Hi. I want to solve the following equation using ode45 solver.
How to set the initial conditions?
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)