Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

1D 2nd-Order nonlinear differential eqn using Matlab

  1. Sep 2, 2005 #1
    Hi all,

    I am more into physics than maths and I need you guys to please help me out. Also I am new to programming in MATLAB. I have done some but still I consider myself a novice. I need to solve the problem in MATLAB.

    So here is the problem.Its a 2 point BVP(boundary value problem) along a straight line ,say along x-axis , from -L to +L.


    A*N''=B+(C*N)+(D*N^2)+f(N)


    [prime denotes differentiation w.r.t x, N~N(x) , f(N) is a simple function of N, say log(N), All constants are known]

    Now I need to solve this BVP for the following boundary conditions:N(-L)=B1 and N(L)=B2, Both B1 and B2 are also known.

    The problem is that this particular equation needs to be solved using the "Relaxation Method" and Runge-Kutta/ Shooting method's are not considered suitable traditionally for this diffusion like problem.

    Can any body suggest me a code for this problem in Matlab. I haven't seen any function in MATLAB capable of solving NL BVP.

    Please help.

    Sincerely,

    JAM
     
  2. jcsd
  3. Sep 2, 2005 #2

    lurflurf

    User Avatar
    Homework Helper

    So here is the problem.Its a 2 point BVP(boundary value problem) along a straight line ,say along x-axis , from -L to +L.


    A*N''=B+(C*N)+(D*N^2)+f(N)


    [prime denotes differentiation w.r.t x, N~N(x) , f(N) is a simple function of N, say log(N), All constants are known]

    Now I need to solve this BVP for the following boundary conditions:N(-L)=B1 and N(L)=B2, Both B1 and B2 are also known.
    First thing I do not see the point of that particular form. Why not consider
    N''=f(N)

    A few things to consider when solving these types of problems
    -How much accuracy is needed?
    what is more important accuracy in the DE of in the BC?
    -How many calculations are you willing to perform?
    Would 1 minute, 1 hr, 1 day ,1 week running on a computer be best?
    -Is any thing known that would simplify/worsen matters?
    Does f(N) have any known properties
    Might there be a boundry layer?
    Do you expect good convergence?

    The problem is that this particular equation needs to be solved using the "Relaxation Method" and Runge-Kutta/ Shooting method's are not considered suitable traditionally for this diffusion like problem.

    Do you mean you wish to use a relaxation method and definately not a shooting method? What relaxation method do you exactly mean?
    I think the idea you have in mind is
    N''=f(N) N(-L)=b1 N(L)=b2
    so you guess N and call your guess N(t,0)
    N should then satisfy a difference equation approximatly, that approximate the DE such as
    N(t+h)-2N(t)+N(x-h)=(h^2)f(N(x)) h a small number
    so you would take your guess and repeatedly improve it by itteration
    solve
    N(t+h,0)-2N(t,1)+N(x-h,1)=(h^2)f(N(x),1)
    for N(x,1)
    if f makes solving hard you could iterate for it or solve instead
    N(t+h,n)-2N(t,n+1)+N(x-h,n+1)=(h^2)f(N(x),n)
    you could also use
    N(t+h,n)-2N(t,n+1)+N(x-h,n)=(h^2)f(N(x),n)
    but results would likely be worse
    or you could change direction
    ie switch between solving
    N(t+h,n)-2N(t,1)+N(x-h,n+1)=(h^2)f(N(x),n+1)
    and
    N(t+h,n+1)-2N(t,n+1)+N(x-h,n)=(h^2)f(N(x),n+1)
    or even use
    N(t+h,*)-2N(t,n+1)+N(x-h,*)=(h^2)f(N(x),1)
    where you update N(t,n) using its two near niehbors, but when which point is updated is not sequential, but some other criteria.

    note: that due to BC we have
    N(-L,n)=b1 N(L,n)=b2 for all n
     
  4. Sep 3, 2005 #3
    problem details

    I am attaching the details of my problem with this reply.
    If you dont get the details can you please email me at jayanta.mukherjee@gmail,com so that I can email you the details of my problem.

    You are absolutely spot on with your suggestions.

    I would like to discuss the problem with you in details.

    thanks in advance.

    jam
     

    Attached Files:

  5. Sep 4, 2005 #4

    lurflurf

    User Avatar
    Homework Helper

    First the logic behind the recipe.
    Boundry value problems are tricky because the differential equation determines the local behavior of the function, while the boundry conditions are at different points. Thus several techniques use the idea of a sequence of functions which satisfy the boundry conditions, that (hopefully) converge to also satisfy the differential equation.
    So we start out with some function that satisfys the BC, if we don't know any thing better we may try a linear function
    N(x,0)=N(-L)+N(L)+(N(L)-N(-L))x/(2L)
    now if the differential equation is
    N''=f(N)
    N should for small h approximately satisfy
    N(x-h)-2N(x)+N(x+h)=(h^2)f(N(x))
    I will assume this is readily solved for N(x)
    N(x)=g(N(x-h),N(x+h))
    if this is not the case the N(x) can be iterated for, or the previous N(x) can be used inside f, in both cases it is hoped convergence still occurs.
    So the idea is we can repeatedly go through our points and improve our guess.
    in psuedo code it might look something like
    define vector N={...};our first guess
    while(dN'dN<eps)
    for i=1,...,n-1
    dN(i)=g(N(i-1),N(i+1))-N(i)
    N(i)+=w*dN(i)
    end(for)
    end(while)
    where w is the damping
    0<w<1 for under relaxation
    w=1 for relaxation
    1<w<2 for over relaxation

    in matlab (I am rusty so this may be off)

    n=128;
    w=1.0;
    loop=0;
    maxloop=100;
    N(1:n)=b1+b2+2*(b2-b1)([0:n]-n/2)/n;
    dN(1:n-1)=1;
    while(loop<maxloop & dN'*dN<eps)
    dN(1:n-1)=\*some function or loop*/;
    N(1:n-1)=N(1:n-1)+w*dN(1:n-1)=;
    loop=loop+1;
    end
    plot(N)
    N

    remember there are many variations
    -different orders of update
    -changing w during looping
    -interpolating new values of N
    -changing the spacing during iteration ie first n=32 then 64 then 128 ect

    I am a bit confused why your boundry conditions seem to depend on N
    k.exp(-N(x)/w)

    The iteration number is just how many times we have updated our guess. The first guess can be k=0 the next k=1 then k=2 ect.

    a few matlab tip pages
    http://www.cyclismo.org/tutorial/matlab/control.html
    http://www.eng.auburn.edu/~sjreeves/Classes/DSP/DSP.html
    http://www.me.pdx.edu/~gerry/MATLAB/programming/performance.html
     
    Last edited: Sep 4, 2005
  6. Sep 4, 2005 #5
    nonliear BVP

    Thanks a million for the help. I am beginning to understand what I need to do. But I still have some questions. Please help me with these.
    Firstly kindly explain, your choice of linear function.
    What is the argument '0' specify in the LHS of the above function?Did you use it to denote the initial N value?
    Secondly, please explain what do mean by these two statements. I don't quite understand the symbols in it.
    Is it the same convergence criteria which I gave in my attachment.If it is, then probably I don't understand the symbols that you have used. What does 'd' stand for?
    Also I didn't understand the statement below. Can you please elaborate. What is it used for?
    As for the boundary condition , I acually have a Neumann type bounday condition at the boundaries.If written in Dirichlet format it should be an exponentially decreasing function of N(I hope I am correct), which is what I want physically as well.The value of N should decrease to zero from its value at the boundary, away from it(on both sides).

    By the way how and where are the boundary conditions taken care off in your code.Kindly explain?

    The confusion about 'k' as an index in my previous post(attachment) is still there. I don't understand why 'k' is with 'i+1' and 'i-1' , whereas 'k+1' is with 'i'?

    Is is possible to derive out a finite difference relation with both 'k' and 'i' in it the way it is given in my previous post.That will help me understand the whole idea better.

    Also your logic behind the statement for the matlab code is confusing for me.
    Can you please explain your idea behind it.

    Thanks for your valuable time.

    Waiting for your reply eagerly.

    Sincerely,

    JAM
     
  7. Sep 5, 2005 #6

    lurflurf

    User Avatar
    Homework Helper

    First your Neumann BC are most unfortunate. They make can make stability and uniqueness a problem. Often Neumann BC do not determine a unique solution and more conditions are needed. Though I could see were no flux BC are appropriate for a diffusion problem.
    I take this to mean
    N'(-L)=0=N'(L)
    or possibly
    N'(-L)=A1 N'(L)=A2

    Thanks a million for the help. I am beginning to understand what I need to do. But I still have some questions. Please help me with these.
    Firstly kindly explain, your choice of linear function.

    Quote:
    N(x,0)=N(-L)+N(L)+(N(L)-N(-L))x/(2L)
    This was when I thought there was hope for Dirichlet BC. In that case and barring a better guess a linear function is a simple function that meets the Dirichlet BC. I guess this can now be ignored.

    What is the argument '0' specify in the LHS of the above function?Did you use it to denote the initial N value?
    The extra argument (0 here) I use as an iteration number. N(x,0) is the first guess. N(x,1) is a refinment of that and so on.
    It is hard to guess the solution well, so just pick something simple that meets the BC and hope for the best. Now that we are dealing with Neumann BC something like
    N(x,0)=0 or
    N(x,0)=1+(N'(L)+N'(-L))x/(2L)+(N'(L)-N'(-L))x^2/(4L)
    should we known flux instead of no flux.
    Again the idea is just to guess a simple function that satisfys the BC

    Secondly, please explain what do mean by these two statements. I don't quite understand the symbols in it.

    Quote:
    while(dN'dN<eps)
    Is it the same convergence criteria which I gave in my attachment.If it is, then probably I don't understand the symbols that you have used. What does 'd' stand for?

    yes it is the same
    dN(x,n)=N(x,n)-N(x,n-1)
    it is the change in N between iterations
    dN'dN is the hermetian inner product
    presuming we mean the discreat version and N real
    dN'dN=dot(dN,dN)=dN^2=sum([N(x,n)-N(x,n-1)]^2;as x takes each chosen value)
    Thus when dN'dN is small each iteration only changes N slightly and we stop
    This can be dangerous if it is a slowly converging sequence, we may be far from the answer. d just stands for difference or change in. That notation is well suited to matlab which understands dN'*dN correctly.

    Also I didn't understand the statement below. Can you please elaborate. What is it used for?

    Quote:
    dN(i)=g(N(i-1),N(i+1))-N(i)

    That just means I assume a solution to
    N(x-h,k)-2N(x,k+1)+N(x+h,k)=(h^2)f(N(x,k+1))
    N(x,k+1)=g(N(i-1),N(i+1))
    If this is hard to solve some aproximation or iteration can be used
    then
    dN(i)=g(N(i-1),N(i+1))-N(i)
    just means
    I get what should be added to N(i) in the itteration
    It is helpful to write
    dN=g(N(i-1),N(i+1))-N(i)
    N=N+dN
    in the loop because
    -it keeps only one N to avoid confusion
    -we need dN anyway for the stoping criteria
    -N=N+w*dN is easily done if it is later decided over relaxation would help

    As for the boundary condition , I acually have a Neumann type bounday condition at the boundaries.If written in Dirichlet format it should be an exponentially decreasing function of N(I hope I am correct), which is what I want physically as well.The value of N should decrease to zero from its value at the boundary, away from it(on both sides).
    The BC is the hard part. I do not know how they would be best delt with. At first guess I would expect having the descrete version satisfy the finite difference equivalent would be better than subing in Dirichlet BC, but I am not sure.

    By the way how and where are the boundary conditions taken care off in your code.Kindly explain?
    As I mentioned before the nice thing about relaxation is in each loop the values returned satisfy the BC. We iterate to try and satisfy the DE. This is a nice change from methods were the DE is easily satisfied, but then the BC are difficult to satisfy. In that sample code I assumed Dirichlet BC, my first guess meet them, and the iterations never changed the boundry values.

    The confusion about 'k' as an index in my previous post(attachment) is still there. I don't understand why 'k' is with 'i+1' and 'i-1' , whereas 'k+1' is with 'i'?
    we discreatize the DE to
    N(i-1)-2N(i)+N(i+1)=h^2f(N(i))
    we make a first guess N(x,0) where x takes the needed values
    we want to satisfy the DE so we keep refining N(x) so that the new value of our point and the old values of its neighbors satisfy the difference equation. With the hope that eventually the BC and approximating Difference Equation are both satisfied (at least approxiamately). so
    given N(i-1,k) and N(i+1,k) (old guesses)
    we want a new guess N(i,k+1) so that N(i-1,k), N(i,k+1),and N(i+1,k) satisfy
    N(i-1,k)-2N(i,k+1)+N(i+1,k)=h^2f(N(i,k+1))

    Is is possible to derive out a finite difference relation with both 'k' and 'i' in it the way it is given in my previous post.That will help me understand the whole idea better.
    Again as above two different things are happening. We have a difference equation and iteration. We pick one point N(i,k) and say that does not satisfy the difference equation, no problem I will just pick a new value N(i,k+1) that does. Then N(i+1,k) does not do so we pick N(i+1,k+1) and so on. The iteration numbers can be confusion because we chang all the points over and over (sometimes in different orders). It is probably better just count the loops.


    Also your logic behind the statement for the matlab code is confusing for me.

    Quote:
    N(1:n)=b1+b2+2*(b2-b1)([0:n]-n/2)/n;

    Can you please explain your idea behind it.
    That was just a way (in matlab) to initalize the vector for the linear equation guess. I did have a typo should be
    N(0:n)=b1+b2+2*(b2-b1)([0:n]-n/2)/n;
    also
    N=b1+b2+2*(b2-b1)([0:n]-n/2)/n;
    might work I forget
    I forget some matlab stuff I think an index 0,1,...,n-2,n is okay, but I am not sure

    Now that we are discussing Neumann BC a few notes
    we can derive difference equationsfor the BC
    3N(-L)=4N(-L+h)-N(-L+2h)-2h*N'(-L)
    3N(L)=4N(L-h)-N(L-2h)+2h*N'(L)
    and apply these each time we reach a boundry.
    Now new psuedo code is in order I think iterating alternating foward and back might be good idea with Neumann BC

    define vector N={...};our first guess

    while(dN'dN<eps)

    for i=1,...,n-1
    dN(i)=g(N(i-1),N(i+1))-N(i)
    N(i)+=w*dN(i)
    end(for)

    dN(n)=(4*N(n-1)-N(n-2)+2*h*A1)/3-N(n)
    n(n)+=dN(n)

    for i=1,...,n-1
    dN(n-i)=g(N(n-i+1),N(n-i-1))-N(n-i)
    N(i)+=w*dN(i)
    end(for)

    dN(0)=(4*N(1)-N(2)-2*h*A2)/3-N(0)
    N(0)+=dN(0)

    end(while)

    (second for loop lift right out if only forward iteration is desired)
    I think you can vectorize reverse iteration in matlab, but you better check (if you desire reverse iteration) or write it as forward iteration.
    It is important to vectorize when possible as matlab is not a compiled lannguage and using unvectorized calculations is dreadful slow. By vectorized I mean
    S=sin([0:pi/10000:2*pi]);
    is much faster than using a for loop
     
    Last edited: Sep 5, 2005
  8. Sep 6, 2005 #7
    Thanks for making things clear for me. I am working on your suggestions. I will come back soon.

    Regards,

    jam
     
  9. Sep 12, 2005 #8
    Hi !

    Thanks for all the help ,Lurflurf. I have been able to solve the problem. But there is still one thing which needs to be sorted out.It is the convergence criteria.

    My while loop in MATLAB is of the form:

    while k<maxiter

    statement..;
    statement..;
    end

    where 'k' is the iteration index

    But the actual convergence criteria which i am not able to implement is

    e=dN'*dN=N(:,k+1)-N(:,k)<0.0001(say);

    Since dN is calculated inside the while loop , I am not able to specify the criteria e<0.0001 against the While statement outside the loop.

    So I have tried the 'if break' statement inside the 'while loop' , but it too doesn't work.

    Also the dN is getting calculated out to be zero.

    pls tell me if you wish to see my code. I hope you will be able to help me out.

    Sincerely,

    JAM
     
  10. Sep 12, 2005 #9
    Hi, I have been able to resolve the matter using the if break staement inside the while loop.
    However I still have a question. Please give me your opinion as to which one is a more correct or appropriate convergence criteria:

    dN'*dN=eps<1e-6(say);
    or
    sqrt(norm(dN))<1e-6(say);

    I will be waiting for your reply.

    Thanks for all the help.

    Jam
     
  11. Sep 12, 2005 #10

    lurflurf

    User Avatar
    Homework Helper

    I think you switched things
    assuming the L2-norm
    sqrt(dN'*dN)=norm(dN)
    so either can be used
    I am confused by this
    e=dN'*dN=N(:,k+1)-N(:,k)<0.0001(say);
    Are you keeping all the iterations?
    If so watch out for ecessive memory use.
    You should be able to use a while with a variable that changes inside the loop, that is the whole point of using a while instead of a for..
     
  12. Sep 13, 2005 #11
    thanks for the reply.
    there was small error while writing e. it is actually,
    e=dN'*dN=N(p,k+1)-N(p,k)<0.0001(say);
    where 'p' denotes the x-space discritization index. I hope its alright now.
    However can you please elaborate on

    "Are you keeping all the iterations?
    If so watch out for ecessive memory use. "

    Thanks again.
    Jam
     
  13. Sep 13, 2005 #12

    lurflurf

    User Avatar
    Homework Helper

    I just wondered how many previous values of N you keep. The usual thing would be to keep one set of N directly and one set implicitly in dN. Sometimes to track how the iterations progress additional N values are kept. Keeping too many previous values can cause memory use problems.
     
  14. Sep 13, 2005 #13

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    So you got this:

    [tex]aN^{''}=b+cN+dN^2+f(N);\quad N(-L)=B_1,\quad N(L)=B_2[/tex]

    You mind supplying the values of f(N), the constants, and the values at the end points? I'd like to look at it in Mathematica or else I'll just wip up some values and see what I come up with. :smile:
     
    Last edited: Sep 13, 2005
  15. Sep 14, 2005 #14
    Actually, you can add another term on the RHS of the form e*N^3.I did not include this earlier as it made matters more complicated for me. 'e' is another constant(very small as compared to b,c and d).

    You can assume B1=B2=0 and let L=100. f(N) can be just cosnt*(N-No), where No is another constant.

    Please do tell me the method you use for the solution. I have used relaxation method. But I have not included the N^3 term yet. I wish to include it eventually.I have used Matlab.I dont know Mathematica.
     
  16. Sep 14, 2005 #15

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Hey Jam, what specific values are you using for a, b, c, d, No, e? Really, I'll try the shooting method first. Is there some reason this method is not applicable?
     
  17. Sep 15, 2005 #16
    Please use the following values:
    L=100, a=3e9 , b=9.375e+15 c=0.2e9, d=1.4e2 , e=12 , No=1e6 .

    The shooting method and the RK4 has been reported to give inaccurate results when the term 'e' is a highly oscillating function(not a constant as in the posed problem).

    In my actual problem 'e' is a highly oscillating function of 'x'. I posed it as a constant just to make matters simple. You can use e=e*f(x). where f(x) is a highly oscillating symmetric function of x, if u want to go the whole way.

    Although shooting method should be fairly accurate I don't know why it gives errors for as 'x' tends to L. i.e near the boundaries.

    Please do tell me the recipe you will use to solve the given nonlinear DE.

    Thanks for taking interest.

    Jam
     
  18. Sep 15, 2005 #17

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Jam . . . I still don't have it dude. What does 3e9 mean? Is that:

    [tex]3e^9[/tex]

    or:

    [tex]3(12)^9[/tex]

    I suspect it's the former and will run with that. :smile: Even though you said the shooting method is unacceptable I'll still try it and wish to determine why and then will pursue yours and Lurflurf's avenue. It's all in the journey. :smile:
     
  19. Sep 15, 2005 #18
    hey, 3e9 is just 3*10^9.(^ denotes to the power) Hope you get it now. Same for others.
    Jam
     
  20. Sep 15, 2005 #19

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Alright. Thanks Jam. So this is what I'll look at:

    [tex](3\cdot 10^9)y^{''}=(9.375\cdot 10^{15})+(0.2\cdot 10^9)y+(1.4\cdot 10^2) y^2+12y^3+(y-10^6);\quad y(-100)=0,\quad y(100)=0[/tex]

    Please correct if not Ok.
     
  21. Sep 15, 2005 #20

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Alright, it's not happening for me via shooting method. But that's Ok Jam. I'll look into Relaxation Method as this is new for me and I'm here for the ride and will rely on a comment I made as one of my first posting in PF:

    When confronted with a problem you can't solve, put it up and work on a simpler related one. So I'll first look at:

    [tex]y^{''}=-\frac{1}{(1+y)^2};\quad y(0)=y(1)=0[/tex]

    via Relaxation Method. Then maybe another one, maybe another, then I'll return to yours. May take a while and anyway I think Lurflurf helped you ok. :smile:
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: 1D 2nd-Order nonlinear differential eqn using Matlab
Loading...