second order nonlinear ode with polynomial terms

Hi everyone! I would resolve the following nonlinear differential equation:

f(Y) + b(Y) (Y')^2 + g(Y) Y'' = A

where Y is a function of x, i.e. Y = Y(x), and

 f(Y) = a1 + a2*Y + a3*Y^2 + a4*Y^3
 g(Y) = b1 f(Y)/Y^3
 b(Y) = c1 (a1 + a2 Y + a3 Y^2)/Y^3

In this example A = cost, but it could be A = A(x). I have no idea how to solve with matlab.. some suggestions? can I use some usual ode-routines?

Thanks in advance for all your support.

Pinco

[Pinco]

up

[Kye Taylor]

What have you tried so far? Do you know you have to write the ODE as a system of first order differential equations?

Also, do you mean that t and x are the same independent variable? Or is x a function of t? Or are there two independent variables? Whoa dude...

[cr]

If A is a function of Y, it needs to be known what function it is. I dont think its possible to generalize solution for any A. Also need to know what is x and is Y' = dY/dx or dY/dt?

[Jan]

@Pinco: Bumping after 2 hours is not useful. When the contributors do not find enough information to create an answer, reading the question again without any additional explanantions, is a waste of time. So please do not bump after 2 hours without showing, what you have done in this time.

[Pinco]

@Kye Taylor @Chandrakanth, : You are right, there was a mistake (I fixed it): A = A(x). In this case, the function A is completely known.

' = d/dx.

I have to find Y(x), thus I have only independent variable x.

Is it then sufficient to solve an ODE system?

@Jan Simon: I apologize for refreshing, but I was very worried about this solution

Answer 1

[Kye Taylor]

You must write the ODE as a system of first-order ODEs. Use the substitution u1 = y and u2 = y'. Then, you'll end up with equations like

u1' = u2

u2' = F(u1,u2)

where F is a function of u1 and u2 (y and y')... Once you have those equations, create a function named F, like

 function uPrime = F(u)

 uPrime(1) = u(2);
 uPrime(2) =  % code for your F(u1,u2)

Note that the input u should be a two-dimensional vector where the comoponent u(1) represents u1, and u(2) represents u2. The output is also a two-dimensional vector, one element for each first-order differential equation in the system above. Such an interface to F is dictated by the requirements of the fsolve function.

[Pinco]

Thanks for this answer. I can use also the classical ode45 or ode23, right? Using fsolve, how can I impose the boundary conditions?

[Kye Taylor]

Ooops, i meant to say ode** ... in fact yeah i'd use ode23s, depending on what you do with those parameters; I've noticed in experimenting with my solution that there's some funny things depending on the parameters and interval of time you solve over.

Good luck... fun system.

[Pinco]

I implemented the function in this way:

 function xprime = myode(t,x)

 global a1 a2 a3 a4 h0 k

 A = (12.* x(1).^3)/h0^2;
 B = a1 + 2*a2.*x(1) + 3*a3.*x(1).^2 + 4*a4.*x(1).^3;
 C = 3*a1 + 4*a2.*x(1).^2 + 3*a3.*x(1).^2;

 xprime = [x(1); ...
        A.*(1- k./B) + (C./x(1).*B).* x(2).^2];

Now I would impose some boundary conditions, i.e. Y(inf) = L, y(-inf) = l, dY(inf)/dx = dY(-inf)/dx = 0. How can I do this?

[Pinco]

Now I'm trying to use bvp4c for boundary value problems. According to the second example in the Matlab help for bvp4c, I wrote this code:

 global a0 a1 a2 a3 a4  h0 k Ja Jb

 a0 = 2.0377638272727268 * 10^3; 
 a1 = -7.105521894545453 * 10^3;
 a2 = 9.234000147272726 * 10^3;
 a3 = -5.302489919999999 * 10^3;
 a4 = 1.1362478399999998 * 10^3;
 h0 = 45.5 * 10^-10;
 k = 5.92; 
 Ja = 1.025;
 Jb = 1.308;

 meshX = linspace(0,16);
 Y0 = [Ja, 0]; %[Y(0), Y'(0)]

 init = bvpinit(meshX,Y0);
  sol = bvp4c(@myode2,@mybc2,init);

======

 function yprime = myode2(x,y)
 global a1 a2 a3 a4 h0 k

 A = (12.* y(1).^3)/h0^2;
 B = a1 + 2*a2.*y(1) + 3*a3.*y(1).^2 + 4*a4.*y(1).^3;
 C = 3*a1 + 4*a2.*y(1).^2 + 3*a3.*y(1).^2;

 yprime = [y(1); ...
          A.*(1- k./B) + (C./y(1).*B).* y(2).^2];

=====

 function res = mybc(ya,yb)
 global Ja Jb

 res =[ya(2)
       yb(2) 
       ya(1) - Ja
       yb(1) - Jb];

I would impose the boundary conditions on the first derivative (equals to 0 at the extremes) and the function itself (equals to Ja and Jb, respectively). Matlab returns some errors.. where is the error?