Interpolation for large number of points

Posted on 2009-06-26 19:51 By manjees (

47)

Categories:

How do I ...? with Maple

I have this code:

ade:= diff(y(r),r,r)+2*diff(y(r),r)/r+486*polylog(3/2,-exp(1-1/43*y(r))) = 0;
ys:= dsolve({de, y(0) = 0, D(y)(0) = 0}, y(r), series, order=10);

yn := convert(evalf(ys), polynom);yn1 := eval(yn, r = .1); dyn1 := eval(diff(yn, r));

dyn2 := eval(dyn1, r = .1);S := dsolve({de, y(.1) = 1.266393114, (D(y))(.1) = 25.18200436}, y(r), numeric);with(plots); odeplot(S, [r, y(r)], r = 0.1e-1 .. 200);

# I get a plot from this step, but I need an expression for y(r) also. A crude method is to get S(1);S(2) and so on and copy paste those values for interpolation.

Question: Would someone mind helping me with a generic code for that step?

Expression

Comment on 2009-06-26 20:14 from Robert Israel ( Maple Rating 4

2434)

There is presumably no closed-form expression for y(r). If you want a good approximation on a particular interval, you might try something from the numapprox package, perhaps chebyshev or chebpade. Interpolation is _not_ likely to give good results. But I don't know why you think you need "an expression" for y(r).

Hmmm, it looks like chebpade has trouble with dsolve(..., numeric) output. This may warrant further investigation.

Why I need expression for y(r)

Comment on 2009-06-26 20:15 from manjees (

47)

I have to plug that expression for y(r) in two simultaneous integral equations and solve them using NRapson. I need to keep changing everything till Nrapson starts converging.

Thanks again for all the help
MS

integral equation

Comment on 2009-06-28 23:22 from Robert Israel ( Maple Rating 4

2434)

Well then, you're asking the wrong question. You don't need an expression for y(r), you need to know how to calculate integrals involving y(r). For example, suppose you want to calculate int(y(r)*cos(r), t = 1 .. 2) where y(r) is a solution of differential equation de with initial conditions ic:

> S := dsolve({de, ic}, y(r), numeric, output=listprocedure);
Y:= subs(S, y(r));

Now the following might work:

> evalf(Int(Y(r)^2 * cos(r), r = 1 .. 2));

If not, you might try this:

> S2:= dsolve({diff(J(r),r) = Y(r)^2 * cos(r), J(1) = 0}, J(r), numeric);
subs(S2(2), J(r));

Both methods work

Comment on 2009-06-29 14:07 from PatrickT ( Maple Rating 4

530)

Both methods work ... but yield slightly different values. Interesting.

> de := diff(y(r), r) = 1:
   ic := y(0)=1:
> S := dsolve({de, ic}, y(r), numeric, output=listprocedure):
   Y:= subs(S, y(r)):

> evalf(Int(Y(r)^2 * cos(r), r = 1 .. 2));

                            0.02404977541

> S2:= dsolve({diff(J(r),r) = Y(r)^2 * cos(r), J(1) = 0}, J(r), numeric):
   subs(S2(2), J(r)); 

                        0.0240498760230187549

not able to get chebyshev workin...

Comment on 2009-06-27 00:16 from manjees (

47)

basically m not able to direct the dsolve output to this command..

any suggestions
thanks

not sure

Comment on 2009-06-27 09:09 from PatrickT ( Maple Rating 4

530)

I'm still not sure I understand your question, but perhaps this thread will help.

from a plot structure named p, you can extract the data points with

map(x->op(1,x),op([1,1],p)):
map(x->op(2,x),op([1,1],p)):

the first line extracts the x-values, the second line the y-values, you can then rearrange them into an array. There is probably a way to do that with one line of code, but since I just found out about this map trick two minutes ago I'm not sufficiently competent to help with this.

All the best manjees.