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

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

lurflurf
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

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 [email protected],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: 222
lurflurf
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:
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

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

Regards,

jam

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

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

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

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

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

saltydog
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:
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.

saltydog
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?

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

saltydog
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. hey, 3e9 is just 3*10^9.(^ denotes to the power) Hope you get it now. Same for others.
Jam

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

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

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

saltydog
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[] = 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[] = 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[] = 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

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

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