I need help solving this second order differential equation?

I was given a second order function that models longitudinal displacements in a longitudinally loaded elastic bar.

a(x)*u''(x) + a'(x)*u'(x)=f(x), 0 <= x <=1

I know that the left end is at x=0 and the right end is at x=1. So that means:

u(0)=u(1)=0 (does not bend at the ends)

I know a(x) represents both the elastic properties and the cross-sectional area of the bar and u(x) is the displacement at point x.

for the first part I am told that a(x)=1+x and f(x)= 5*sin(2*pi*x)^2 (this is the case where the force is applied symetrically and is strongest at x=.25 and x=.75)

I need to use dsolve to solve for the problem and plot the result on the interval [0,1]

This is what I tried to do:

ode1= 'D2y*(1+x)+ Dy= 5*sin(2*pi*x)^2';

sol=dsolve(ode1,'y(0)=0','y(1)=0');

ezplot(sol)

Howver the graph I got from this just does not look right. Can someone please help me out?

[Robin]

...And oh yea, the graph i get is a x-t graph for some reason??

[bym]

t is the default independent variable, try using:

sol=dsolve(ode1,'y(0)=0','y(1)=0','x');

you'll get a nasty expression involving exponential integrals that EZ plot can't handle

[Robin]

I actually did try that, and that's right I got a long nasty expression, but the thing is there are some 'Ei' and 'i' in there. So are those 'Ei' just constants? What if i need to plot it? how would I go about plotting it?

[Walter Roberson]

Ei would be the Elliptic Integral. It is a constant if all of the parameters to it are constants.

i would be the square root of negative 1.

[Robin]

ok yea I figured the i as a complex thing but I am not at all familiar with "elliptic integral"...so how would I deal with this? All I need to do is plot the solution on an interval from 0<x<1, so I was thinking I could just make an array of X and Y points, and just use plot(X,Y). But I am just not sure about what to do about the "elliptic integral" part?

[Robin]

ok yea I figured the i as a complex thing but I am not at all familiar with "elliptic integral"...so how would I deal with this? All I need to do is plot the solution on an interval from 0<x<1, so I was thinking I could just make an array of X and Y points, and just use plot(X,Y). But I am just not sure about what to do about the "elliptic integral" part?

[bym]

Ei is short for exponential integral and i is sqrt(-1)

[Walter Roberson]

http://en.wikipedia.org/wiki/Elliptic_integral

The expression can be converted to the equivalent hypergeometric form,

y(x) = (-5*(-ln(8*Pi)+ln(4*Pi))*(1+x)*hypergeom([1/2],[3/2, 3/2],-4*(1+x)^2*Pi

^2)+(-10*ln(4*(1+x)*Pi)+10*ln(4*Pi))*hypergeom([1/2],[3/2, 3/2],-16*Pi^2)+(5*

ln(4*(1+x)*Pi)-5*ln(8*Pi))*hypergeom([1/2],[3/2, 3/2],-4*Pi^2)+5*ln(4*(1+x)*Pi

)+(-5+5*x)*ln(4*Pi)-5*x*ln(8*Pi))/(-2*ln(8*Pi)+2*ln(4*Pi))

[bym]

you might be able to plot is using subs. You can evaluate the exponential integral using the following:

mfun('Ei',-4i*pi)

ans =

-0.0061 - 3.0630i

[Robin]

ok so I tried to replace what I had for that long solution with mfun, but I get an error:

>> y=(5*x)/2 - (5*log(x + 1))/2 + (5*i*mfun('Ei',((-4)*pi*i)) - 5*i*mfun('Ei',(4*pi*i)))/(16*pi) - (5*i*mfun('Ei',((-4)*pi*i*(x + 1))))/(16*pi) + (5*i*mfun('Ei',(4*pi*i*(x + 1))))/(16*pi) - (log(x + 1)*(40*pi - 40*pi*log(2) + 5*i*mfun('Ei',((-4)*pi*i)) - 5*i*mfun('Ei',(4*pi*i)) - 5*i*mfun('Ei',((-8)*pi*i)) + 5*i*mfun('Ei',(8*pi*i))))/(16*pi*log(2))

??? Error using ==> str2num at 33

Requires string or character array input.

Error in ==> mfun>computeX at 72

x = str2num(x); %#ok

Error in ==> mfun at 43

[x,siz,nans] = computeX(varargin{:});

[Walter Roberson]

My mistake earlier: Ei is the Exponential Integral, not the Elliptic Integral. The hypergeometric conversion I showed is still valid, though.

I do see any immediate reason why you would be getting the str2num errors, but I can make the suggestion that you optimize using

t = mfun('Ei', pi*i*[-4, 4, -8, 8, -4*(x+1), 4*(x+1)]);

and then replace the individual mfun calls with t(1), t(2), t(5, t(6), t(1), t(2), t(3), and t(4) respectively.

Answer 1

[Neels]

My suggestion would be to rearrange the equation as

u''(x) =f(x)/a(x)- u'(x)/a(x)

z = f(x)/a(x)-y/a(x) , where y is first derivative of u(x) and z is first derivative of y. Then solve it using ode45, giving the initial values. I think matlab computes faster using first order derivatives

[Robin]

I did try to rearrange the equation every way I could think of, but I still get the same nasty result with "Ei" as part of the answer. We have to solve it symbolically, using dsolve, the problem specifically says to do that :(