Solve first order differential equation boundary value problem using Matlab

Question

I am trying to solve a model of a physical structure with the following equation:

$y'(x) + \frac{-ax+b}{-x^2+x+m}y(x) + \frac{k}{-x^2+x+m}=0$

The boudnary condition is:

y(0) = -y(1)

where $a = 160$, $b = 6500$, $m = 30$ and $k = 700$ are positive constants. I am only interested in [0, 1].

I have tried to put it into Matlab Mupad or Wolframalpha; however, the solution could not be described as elementary mathematics (it showed something as hypergeometric function). I would like to look at the solution numerically. Since this is not an initial value problem, I do not think ode45 is a good solver in this case. I have googled bvp4c - boundary value problem solver of Matlab. Unfortunately, all of them are about two-point second order ODE. I wonder if someone can give me a hint or guidance how to do it. The preferred software is Matlab, but I am fine with other software also.

Answer 1

Since the ODE was changed after my preceeding answer, I open a second answer, but I don't suppress the first one because it is anyways a interesting exercise.

I have not MATLAB at hand but any other soft for numerical solving of differential equation can be used with the method below.

One can see that $y'=0$ when $(-ax+b)y+k=0$ , hense $y=-\frac{k}{-ax+b}$ . With the example of data and for $0<x<1$, the values of $y(x)$ for $y'=0$ are in the range $-\frac{700}{-160+6500} < y < -\frac{700}{6500}$.

If we start with a value of $y(0)$ not too far from these values, after $y(x)$ comes at the point where $y'=0$ then $y(x)$ remains constant.

So, the process is very simple : Set $y(0)= -\frac{700}{6500}$ for example (other values not too far would lead to almost the same result). The numerical calculus gives $y(1)=-0.110397$ . This is already close to $-y(0)$

In order to obtain a better accuracy, restart the numerical calculus with $y(0)=0.110397$ . The result is $y(1)=-0.11039691$ . This is the wanted result $y(1)\simeq -y(0)$.

The graphical representation of the calculus of $y(x)$ is shown on the next figure :

Only the part of the curve for low values of $x$ is represented since $y(x)$ remains constant up to $x=1$.

An elementary algorithm :

The time of computation and drawing the curve $y(x)$ is less than 1 second on an ordinary PC.

Answer 2

Analytic solution : $$y(x)=-\int \frac{-ax+b+k}{-x^2+x+m} = \frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m | +C$$

With condition $y(0)=-y(1)$ :

$y(0)=\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|m | +C$

$C=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|$

$$y(x)=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|+ +{\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m |} $$

As a consequence :

$y(1)=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m| + {\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|m|} $

$y(1)=y(0)+2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$

The second condition $y(1)=-y(0)$ implies :

$-y(0)=y(0)+2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$

$y(0)=-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$

$$y(x)=-2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|+ +{\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m |} $$

Answer 3

Here is some Matlab code to solve this problem. It's basically a shooting method that uses root-finding to find the appropriate initial condition. There is no guarantee that such an initial condition exists, and if it does, whether or not it would be unique.

The first part of the code defines the ODE you want to solve:

a=160
b=6500
m=30
k=700
der=@(x,y) ((-a*x+b).*y+k)./(x.^2-x-m);

Now set up some code that, when given a value for $y(0)$, solves the ODE and returns $y(0)$ and $y(1)$.

f=@(T) T.y([1 end]);
sol=@(y0) f(ode45(der,[0 1],y0));

Use the Matlab builtin fzero function to find when $y(0)+y(1)=0$,

Y0=fzero(@(y0) sum(sol(y0)),0.5)

and use that $y(0)$ to solve the ODE properly and plot the solution:

[x,y]=ode45(der,[0 1],Y0);
plot(x,y)

This gave me $y(0)=0.11037$ and $y(-1)=-0.11037$.