# 1D 2nd-Order nonlinear differential eqn using Matlab

• MATLAB
jam_27
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.

Sincerely,

JAM

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

jam_27
problem details

lurflurf said:
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

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.

jam

#### Attachments

• My equation is.doc
47.5 KB · Views: 322
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:
jam_27
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.
N(x,0)=N(-L)+N(L)+(N(L)-N(-L))x/(2L)
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.
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?
Also I didn't understand the statement below. Can you please elaborate. What is it used for?
dN(i)=g(N(i-1),N(i+1))-N(i)
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.
N(1:n)=b1+b2+2*(b2-b1)([0:n]-n/2)/n;

Sincerely,

JAM

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
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;

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:
jam_27
Thanks for making things clear for me. I am working on your suggestions. I will come back soon.

Regards,

jam

jam_27
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

jam_27
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);

Thanks for all the help.

Jam

Homework Helper
jam_27 said:
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);

Thanks for all the help.

Jam
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..

jam_27
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

Homework Helper
jam_27 said:
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
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.

Homework Helper
So you got this:

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

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.

Last edited:
jam_27
saltydog said:
So you got this:

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

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.

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.

Homework Helper
jam_27 said:
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.

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?

jam_27
saltydog said:
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?

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

Homework Helper
jam_27 said:
L=100, a=3e9 , b=9.375e+15 c=0.2e9, d=1.4e2 , e=12 , No=1e6 .

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

$$3e^9$$

or:

$$3(12)^9$$

I suspect it's the former and will run with that. 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.

jam_27
hey, 3e9 is just 3*10^9.(^ denotes to the power) Hope you get it now. Same for others.
Jam

Homework Helper
jam_27 said:
hey, 3e9 is just 3*10^9.(^ denotes to the power) Hope you get it now. Same for others.
Jam

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

$$(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$$

Homework Helper
saltydog said:
Alright. Thanks Jam. So this is what I'll look at:

$$(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$$

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:

$$y^{''}=-\frac{1}{(1+y)^2};\quad y(0)=y(1)=0$$

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.

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

$$(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$$

Thats perfect.
Jam

jam_27
saltydog said:
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:

$$y^{''}=-\frac{1}{(1+y)^2};\quad y(0)=y(1)=0$$

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.

This particular problem (the one which you wish to try before going into mine) is described in the following link

http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html

I like the way this guy described the relaxation method. But got lost a bit while trying to understand the recipe. You see I am not good at mathematical jugglery.

If you get hold of the recipe, Please do elaborate it to me.
Meanwhile I worked according to Lurflurf's suggestions and got a 'hopefully' correct result.

Jam

Homework Helper
You know what Jam, I just find this amazing. So using the reference you gave regarding Relaxation Method, I wrote the Mathematica code below to solve the BVP:

$$y^{''}=-\frac{1}{(1+y)^2};\quad y(0)=y(1)=0$$

What I find amazing is that I can use just random numbers as the initial guess! The first plot exhibits my initial random guess. Next I run the code with deltax=0.1 a total of 5 times. The second plot shows how the initial random guess is "relaxed" towards the final solution. The third plot is a shooting-method solution superimposed with the fifth iterate of the Relaxation Method. I'm surprised at how close the results are at just such a coarse resolution of deltax and number of iterates.

Code:
<< LinearAlgebraTridiagonal
(* Initialize variables *)

mtotal = 11;
xstart = 0;
xend = 1.0;
deltax = (xend - xstart)/(mtotal - 1);
y = Table[0, {mtotal}];
e = Table[0, {mtotal}];
a = Table[0, {mtotal}];
b = Table[1, {mtotal}];
c = Table[0, {mtotal}];
plotlist = Table[{0, 0}, {mtotal}];
e[[1]] = 0;
e[[mtotal]] = 0;
rhs[yval_] := -1/(1 + yval)^2;

(* Load y array with initial guess of solution *)

xval = xstart;
For[i = 1, i <= mtotal, i++,
y[[i]] = Random[Real, {-1, 1}];
plotlist[[i]] = {xval, y[[i]]};
xval += deltax;
];
y[[1]] = 0;
y[[mtotal]] = 0;
plotlist[[1, 2]] = 0;
plotlist[[mtotal, 2]] = 0;
(* save plot of initial guess *)
pt1 = ListPlot[plotlist, PlotJoined -> True]

(* Calculate error vector *)

For[i = 2, i <= mtotal - 1, i++,
e[[i]] = y[[i + 1]] - 2 y[[i]] + y[[i - 1]] - deltax^2  rhs[y[[i]]];
];

(* Calculate Taylor expansion of error vector up to linear terms *)

For[i = 2, i <= mtotal - 1, i++,
a[[i]] = 1;
b[[i]] = -2(1 + deltax^2/(1 + y[[i]])^3);
c[[i]] = 1;
];

(* Now set up tridiagonal matrix Ax = -e and solve for delta y(i) *)

adiag = Table[a[[i]], {i, 2, mtotal}];

bdiag = Table[b[[i]], {i, 1, mtotal}];

cdiag = Table[c[[i]], {i, 2, mtotal}];
cdiag[[1]] = 0;

deltayVector = TridiagonalSolve[adiag, bdiag, cdiag, -e]

(* Now update y values *)

For[i = 2, i <= mtotal - 1, i++,
y[[i]] += deltayVector[[i]];
];

(* continue iterating until staisfied with results *)

#### Attachments

• showplot1.JPG
6.6 KB · Views: 467
• showplot2.JPG
8.6 KB · Views: 469
• showplot3.JPG
4.7 KB · Views: 480
jam_27
ya that’s true , the relaxed solution hardly depends on the initial guess. I also find the same behavior in my results.

But I don't know Mathematica at all. I work with Matlab and C++.

Can you please describe me the recipe that you understood from the reference link that I gave in my last post. I am not clear about certain things in that web page. Can you please describe the recipe in your own words? It will be a great help.

Thanks
Jam

Homework Helper
jam_27 said:
ya that’s true , the relaxed solution hardly depends on the initial guess. I also find the same behavior in my results.

But I don't know Mathematica at all. I work with Matlab and C++.

Can you please describe me the recipe that you understood from the reference link that I gave in my last post. I am not clear about certain things in that web page. Can you please describe the recipe in your own words? It will be a great help.

Thanks
Jam

I though since your problem was a special case of an elliptic ODE, the next one I'd work on would be:

$$y^{''}=-y+y^3;\quad y(0)=0,y(8)=0.1$$

with:

$f(y)=-y+y^3[/tex] Here's my summary of it. It helps me understand it too as this is all new to me: Keeping [itex]\Delta x$=0.1, we'll have 80 equations in 80 unknowns for the $\Delta y$ matrix.

So we convert the ODE to a difference equation:

$$y_{i+1}-2y_i+y_{i-1}=(\Delta x)^2 f(y_i)$$

And thus the error term for the i-th value is:

$$e_i=y_{i+1}-2y_i+y_{i-1}-(\Delta x)^2 [-y_i+(y_i)^3]$$

Keeping in mind the error terms at the boundary are 0 and thus much of the looping goes from 2 to maxpoint-1.

Letting $e_i$ be a function of the three y terms:

$$e_i=h(y_{i+1},y_i,y_{i-1})$$

We expand each $e_i$ as a Taylor series up to the linear terms and calculate:

\begin{align*} e_i(\mathbf{Y}+\Delta\mathbf{Y}) &= h(y_{i+1}+\Delta y_{i+1},y_i+\Delta y_i,y_{i-1}+\Delta y_{i-1}) \\ &= e_i+\frac{\partial e_i}{\partial y_{i+1}} \Delta y_{i+1}+ \frac{\partial e_i}{\partial y_{i}} \Delta y_{i}+ \frac{\partial e_i}{\partial y_{i-1}} \Delta y_{i-1} \\ &= e_i+\alpha_i \Delta y_{i-1}+ \beta_i \Delta y_{i}+ \gamma_i \Delta y_{i+1} \\ &= e_i+\Delta y_{i-1}+\left(-2-(\Delta x)^2[-1+3(y_i)^2]\right)\Delta y_i+\Delta y_{i+1} \\ \end{align}

We wish to calculate the $\Delta y_i$ terms such that all the error terms go to zero. Therefore in matrix form we wish:

$$\mathbf{A}\mathbf{Y}+\mathbf{E}=\mathbf{0}$$

or:

$$\mathbf{A}\mathbf{Y}=-\mathbf{E}$$

That's the problem you're having with this right? Took me hours to understand this part. The best way to understand it is to work through the error terms manually, set up a few of the matrix terms manually and see how that matrix equation develops.

If you calculate a few $e_i$ manually, you'd find the matrix A is sparse along the diagonal. In fact it's tridiagonal with the subdiagonal equal to the $\alpha_i$ terms, the diagonal equal to the $\beta_i$ terms, and the superdiagonal equal to the $\gamma_i$ terms.

Now I don't know about you and with some difficulty I could probably code the solution of this matrix equation, Mathmatica just happens to have a tridiagonal solver which I used in the form:

ydeltavalues=TridiagonalSolve[subdigaonal, diagonal, superdiagonal,RHS]

or in terms of the code variables:

The next step is to add these $\Delta y_i$ to each $y_i$ and re-iterate until hopefully the error vector reduces to an acceptable value.

Last edited:
jam_27
I am working on the things you suggested.I am doing the stuff manually.Lets see I get anything. I will come back to you soon.
Thanks,
Jam

Homework Helper
I do wish to follow-up on my investigation into solving:

$$y^{''}=-y+y^3;\quad y(0)=0,\quad y(8)=0$$

via Relaxation method.

The first thing to note is that there are two solutions to this BVP depending if the initial slope is positive or negative. These solutions are shown in the first set of plots (via shooting method). The third plot in the first set is a typical initial random guess (skewed to capture the first solution) to use for the relaxation method. The results are disappointing: Sometimes it converges to the solution as shown in the second set, sometimes it converges to something wrong as shown in the third set. This appears to be the case even if deltax is smaller.

Regarding the problem originally posed in this thread: If I'm unable to obtain reproducible results with this relatively simple problem or even worst, spurious results as shown by the third set of plots, I would have absolutely no faith in the results obtained for more complex equations.

I suspect there are convergence criteria and limitations with this numerical method that a comprehensive study would reveal. However, for the time being, can anyone explain to me why this happens?

Edit: Oh yea, suppose I could be implementing the algorithm incorrectly too.

#### Attachments

• relax3.JPG
16.6 KB · Views: 366
• relax2.JPG
17.2 KB · Views: 418
• relax1.JPG
6.2 KB · Views: 399
Last edited:
Homework Helper
Alright, I'll assume there is some special difficulty with the problem I worked on above but still wish to at least approach the original problem:

$$(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$$

So I'll divide out the coefficient on the derivative and obtain:

$$y^{''}=\frac{9.375}{3}\cdot 10^{6})+\frac{0.2}{3}y+\frac{1.4}{3}\cdot 10^{-7} y^2+4\cdot 10^{-9}y^3+\frac{1}{3}\cdot 10^{-9}y-10^{-3};\quad y(-100)=0,\quad y(-100)=0$$

And in fact for starters, do this one:

$$y^{''}=25+\frac{0.2}{3}y+\frac{1.4}{3}\cdot 10^{-7} y^2+4\cdot 10^{-9}y^3+\frac{1}{3}\cdot 10^{-9}y-10^{-3};\quad y(0)=0,\quad y(22)=0$$

And then gradually ramp up the first coefficient and spread out the interval.

So I run the third one through the relaxation method 3 times and come up with the first plot below. I next fit the data to a 8th degree polynomial and then calculate the derivative at 0. This allows me now to solve the problem as an IVP. The results of that solution, superimposed upon the fitted data, is shown in the second plot. There is general agreement although not of a precise nature.

This is shown in the third plot below which is a plot of the back-substitution (LHS-RHS of the ODE) of the fitted data. The largest difference is about 0.04 which considering the plot goes down to -300 is at least to a coarse approximation ok.

I really think my problem with all this is two-fold: I need a good book on Numerical Methods and I don't have a good handle on precision in Mathematica.

#### Attachments

• relax1b.JPG
5.3 KB · Views: 384
• relax2b.JPG
5.2 KB · Views: 397
• relax3.JPG
3.9 KB · Views: 372
Last edited:
jam_27
Hey Salty, your results seem to be OK but not what I get with my algo. I have verified my result for the originally posed problem, physically, and found the result to match exactly with what is wanted. also my results match the results obtained by commercial FEM software’s which use Newton-Raphson's method to solve this non linear problem. Although I have not included the 3rd order non-linearity, I am getting physically correct results. Can you please tell me the convergence criteria that you are using? May be we can discuss the recipe you are using in your algo and crack it up.

Jam

Homework Helper
jam_27 said:
Hey Salty, your results seem to be OK but not what I get with my algo. I have verified my result for the originally posed problem, physically, and found the result to match exactly with what is wanted. also my results match the results obtained by commercial FEM software’s which use Newton-Raphson's method to solve this non linear problem. Although I have not included the 3rd order non-linearity, I am getting physically correct results. Can you please tell me the convergence criteria that you are using? May be we can discuss the recipe you are using in your algo and crack it up.

Jam

1. What is the physical problem?

2. Can you post your solution (a JPEG plot) to this particular problem?

3. Can you post your solution to the original problem as well as the precise form of the ODE you used and the boundary conditions?

4. Can I access via the net, the commercial FEM software you're using?

I don't have a handle on convergence criteria. I assumed the error vector used in the algorithm would get smaller and smaller as the criteria but I've not seen that reproducibly. I just run the algorithm a few times (for now at least).

jam_27
1. What is the physical problem?

2. Can you post your solution (a JPEG plot) to this particular problem?

3. Can you post your solution to the original problem as well as the precise form of the ODE you used and the boundary conditions?

4. Can I access via the net, the commercial FEM software you're using?

Hey Salty,

I hope this helps:
1) The physical problem is diffusion of electrons in a semiconductor laser.
2) Check attachment
3) Go through my first post , you will find a doc file for details & I have already given you the values of the constants.
4) The software’s are Ansys and FEMLAB. I have used both. Unfortunately they are not for free d/l. You can also use MATLAB's PDETOOL, if you have it.

Thanks for taking interest.

PS: Please do tell me the recipe that you are using. My recipe is given in the doc file. I think you are using a different one.

jam

#### Attachments

• diffusion.jpg
33.3 KB · Views: 427
Homework Helper
Thanks for that plot Jam. I'll try and work towards it.

Yea, your *.doc reference above. That's part of my problem I suppose. I don't have Microsoft Office so I couldn't ever view it. But that's Ok. I'll figure out how to do so one way or another.

Edit: Ok I downloaded a word viewer from the Microsoft site. I'll proceed . . .

The algorithm (recepie) that I use implements the method explained in the reference you gave above.

http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html

I believe I can obtain your result but need to get a better handle on high precision in Mathematica.

Last edited:
Homework Helper
Jam, I do believe I'm making progress with this. I noticed in your doc file that the first coefficient was negative. That changes things. This is the precise equation I solved:

$$y^{''}=-\frac{9.375}{3}\cdot 10^6+\frac{0.2}{3}y+\left(\frac{1.4}{3}\cdot 10^{-7}\right)y^2+\left(\frac{1}{3}\cdot10^{-9}\right)y-\frac{1}{3}\cdot 10^{-3};\quad y(-100)=y(100)=0$$

I've attached a plot of the results. As you can see it's not exactly like your solution. Is there a chance I'm not using the same BVP you used to generate that plot? Please let me know.

Thanks,
Salty

#### Attachments

• squarewave1.JPG
7.4 KB · Views: 381
jam_27
Ok ya, I am sorry about the -ve sign. Its important and I missed it while explaining things to you. Sorry for that. And yes, I forgot to tell you one more thing. The first coefficient (with the -ve sign) on the RHS actually varies from x=-50 to x=50. That’s why I get the plot the way I get it. I am sure you will get the same plot now. By the way your solution looks perfect to me.

PS: Request: can you please explain in your own words the recipe you are using as per the following link. I didn't quite get the hang of it.
http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html

Thanks,
jam

Homework Helper
jam_27 said:
Ok ya, I am sorry about the -ve sign. Its important and I missed it while explaining things to you. Sorry for that. And yes, I forgot to tell you one more thing. The first coefficient (with the -ve sign) on the RHS actually varies from x=-50 to x=50. That’s why I get the plot the way I get it. I am sure you will get the same plot now. By the way your solution looks perfect to me.

PS: Request: can you please explain in your own words the recipe you are using as per the following link. I didn't quite get the hang of it.
http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html

Thanks,
jam

Good Jam. I do believe you've taught me something.

Jam, the method relies on the Taylor expansion of the error terms. Above in post #25, I showed how each error term for each value of y calculated depends on three variables:

$$e_i=F(y_{i+1},y_i,y_{i-1})$$

We wish to calculate the value of this function (for each $e_i$)
for small changes in the variables as per posting #25 above:

\begin{align*}e_i(\mathbf{Y}+\Delta\mathbf{Y}) &= h(y_{i+1}+\Delta y_{i+1},y_i+\Delta y_i,y_{i-1}+\Delta y_{i-1}) \\ &=e_i+\frac{\partial e_i}{\partial y_{i+1}} \Delta y_{i+1}+\frac{\partial e_i}{\partial y_{i}} \Delta y_{i}+\frac{\partial e_i}{\partial y_{i-1}} \Delta y_{i-1} \\ &= e_i+\alpha_i \Delta y_{i-1}+ \beta_i \Delta y_{i}+ \gamma_i \Delta y_{i+1} \\ &= e_i+\Delta y_{i-1}+\left(-2-(\Delta x)^2[-1+3(y_i)^2]\right)\Delta y_i+\Delta y_{i+1} \\\end{align}

We treat the $\Delta y_i's$ as unknowns and set up the n equations so that each is set to zero (that will make each error term zero). Since each term includes an [itex]e_i[/tex] term as per the taylor series of any function, we put those e terms on the right hand side and end up with the matrix equation:

$$\mathbf{A}\mathbf{Y}=-\mathbf{E}$$

The matrix of the coefficients of the delta y terms is tridiagonal because each e'i depends only on its nearest neighbors, the y+1, y, and y-1 variables. That's why the non-zero terms in the matrix go down diagonally. Right?

Solving this tridiagonal matrix equation gives us these values which we then add to the current values of y. Since there is a slight error associated with any taylor expansion, the results are not perfect. That's why iteration is necessary. Hope that helps. If not, which precise part of the algorithm is causing problems for you to understand?

Also, in regards to:

1) The physical problem is diffusion of electrons in a semiconductor laser.

Can you explain how you experimentally verified that the diffusion follows the solution to the BVP? For example, could the output be monitored by an O-scope or some other device which reflects the square-wave phenomenon? Just curious.