4th order ODE, boundary conditions problem

24 views (last 30 days)
I am trying to solve a 4th order linear ODE with boundary conditions. I'm just learning to use bvp4c. The code is the following:
function mat4bvp
solinit = bvpinit(linspace(0,1,5),[1 0 0 0]);
sol = bvp4c(@mat4ode,@mat4bc,solinit);
x = linspace(0,1);
y = deval(sol,x);
plot(x,y(1,:));
function dxdy = mat4ode(x,y);
eps=0.1;
dp=[y(2)
    y(3)    
y(4)    
2.*pi^2.*y(3) +(1./eps).*y(2)-pi^4.*y(1)];
function res = mat4bc(ya,yb);
res=[ya(1)     
yb(1)     
yb(3)     
yb(4)];
But I get the following error message:
??? Error using ==> vertcat?CAT
arguments dimensions are not consistent.
Error in ==> 4odecode>mat4ode at 13
dp=[y(2)
?? Error in ==> bvparguments at 122    
testODE = ode(x1,y1,odeExtras{:});
Error in ==> bvp4c at 120
[n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] =
...
Error in ==> 4odecode at 4?
sol = bvp4c(@mat4ode,@mat4bc,solinit);

Accepted Answer

Andrew Newell
Andrew Newell on 4 Feb 2011
In the assignment to dp (which, by the way, should be dxdy), the fourth line has a space in it which MATLAB thinks is separating two elements. With a couple more corrections, the code is
function mat4bvp
solinit = bvpinit(linspace(0,1,5),[1 0 0 0]);
sol = bvp4c(@mat4ode,@mat4bc,solinit);
x = linspace(0,1);
y = deval(sol,x);
plot(x,y(1,:));
function dxdy = mat4ode(~,y)
eps=0.1;
dxdy=[y(2)
y(3)
y(4)
2.*pi^2.*y(3)+(1./eps).*y(2)-pi^4.*y(1)];
function res = mat4bc(ya,yb)
res=[ya(1)
yb(1)
yb(3)
yb(4)];

More Answers (3)

Gualtiero
Gualtiero on 4 Feb 2011
The solution seem to be extremely sensitive to the form that I give to solinit. Is there any general method to impose a suitable solinit?

Gualtiero
Gualtiero on 4 Feb 2011
in particular, my 4th order ode is
y'''' -2*(pi^2)*y''-(1/eps)*y'+(pi^4)*y=0
b.c.: y(0)=y(1)=0; y''(1)=0; y'''(1)=0
it is not trivial to me to pose an initial guess that would give me the right answer...

Gualtiero
Gualtiero on 4 Feb 2011
edit: an initial guess that satisfy the b.c. is
y=sin(x*pi)*cos(x*pi)+2*(pi^3)*(-(1./3)*(x^3)+(x^2)+(2./3)*x)
However, changing the number of meshpoints in linspace changes drastically the shape of the solution...

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!