my differential equation general solution form requires very high or low values, exceeding double precision capability, while I know its solution is finite

as a brief explanation, I'm trying to solve a first order bessel equation which is part of a more complicated PDE problem, whose solution is as following form:

now we have a bunch of these functions and not one, so there are for example 10 integral coefficients like c1(1), c1(2),... and c2(1),c2(2),....

now when I apply the boundary conditions, a set of linear equations are made by which I can calculate c1, c2,... . in the form "A*C=B".

now as x is near zero, I1(x) goes to zero and K1(x) goes to inf, and as x goes up, I1(x) goes to inf and K1(x) goes to zero. given that y(x) is finite and something less than 0.07, c1 and c2 go to inf or zero depending on their factor. and here the problem begins. working with very high and low values decrease accuracy and might lead to inf or zero.

any idea to overcome this issue?

[J. Alex Lee]

if you can express any equations as ratios of bessel functions, maybe you can switch to approximate expressions at the low and high ends by taking the limits of the ratios of the taylor expansions about 0 and infinity.

[David Goodmanson]

Hi hosein,

how are the boundary conditions determined? Is it a case of making y and x^n dy/dx (for a particular n) to be continuous across boundaries?

[hosein Javan]

hello Mr. Goodmanson. the boundary condition is that the actual 2d function "f" and its derivative must be continous at specific boundaries. but I don't think that is relevant. what matters is that the bessel functions grow very fast and that leads to very large computations that causes unexpected inaccurate results.

you are familiar with this laplace equation, I think what you reffered, I have mentioned in Mr. Lee reply.

in this post I have explained a lot more of the problem in the comments, if you can help. thank you.

https://www.mathworks.com/matlabcentral/answers/580290-how-to-make-symbolic-computations-fast

[hosein Javan]

J. Alex Lee. Mr Lee. thanks for replying. I think you are reffering to a good point that I don't know. I'll ask a question regarding this. consider this example:

we have a differential equation that has this solution:

this solution behaves as the case with bessel function that I mentioned. meaning that r^n and r^-n grow fast. we can "control" this behaviour by multiplying some constants:

where:

note that R_1 and R_2 are constants and since there are "n" individual number of equations, R_1^n and R_2^n can be treated as are constants in every equation.

this way we have ensured that at high "n", the values limits are zeros.

but in this case that we are having bessel functions we cannot do the same:

multiplying by "normalizing" constants:

see that unlike before we have to calculate two large/small numbers then multiplying which results inaccurate. unless as you said there is a limit approximation.

[David Goodmanson]

hello sir,

I still have an inquiry about the derivative aspect. If you need, say, dI/dx and dK/dx, how are those being calculated?

[hosein Javan]

David Goodmanson. hello again sir. I think this is all you need to know. the derivate of I_1 and K_1 can be rewritten in terms of I_0 and K_0.

you're right. I'll give a more thorough explanation. the PDE is an inhomogeneous laplace equation in 2d plane:

where R(x) is a periodic function at x with its main period equal to "L". and at y1<y<y2.

the solution is:

where:

"L" is the (minimum) length of periodicity of f(x,y) which occurs in x coordinate.

"n" is the harmonic number.

the "" is a particular solution.

the "" and "" are integral coefficients and are unknown. these are also y-dependant but in each area they are constants. as following:

by applying boundary conditions we can find the above integral coefficients. (it is quite complicated to explain here but in one sentence it concerns the function and its derivative continuity at each y0 to y6).

now after this we have "n" systems of linear equations in the form:

the elements of A matrix are one of like these:

the goal of the program is to:

1. compute the integral coefficients matrix [C_n]
2. compute the solution using the general solution and the integral coefficients [C_n]

the problem appears in both steps when computing and the rest (the second kind and the zero-order bessel functions). especially at high harmonic numbers of "n", and high values of "p".

-------------------------------

• as for unit scaling, I have already tried that. normally both x and y coordinates are in term of meters. now if we change the unit of y, then the unit of L scales the same as well and the result does not change. therefore the input arguement of the bessel functions is dimensionless. unless it is in terms of radians and maybe we can replace by something? not sure.
• there's another thing I'm thinking that I'm not sure. we wanted the accuracy of f(x,y) at first place and not its bessel functions. maybe the accuracy of bessel functions and their coefficients when they multiply, cancel each other and result into a good accuracy of "f".

-------------------------------

however rigth now I can not set a value lower than L=0.5 when trying to solve for n=99 harmonics.

I hope this time I've explained exactly. sorry again for too much writing.

[J. Alex Lee]

I'm not suggesting "controlling" growth by multiplication with constants...I just off-handedly commented about ratios because, for example, ratios of of n differing by 1 (for same argument v) seem to come up in a lot of fields and so is discussed a lot. Since all explode for large ν, it's hard to get values of their ratios; matlab's built-in besseli seems to fail for . If you only care about ratios, you can come up with their asymptotically good estimates as I mentioned above.

Now in your case above you have and , but for large x, , so maybe you can think about their product. WolframAlpha tells me that the series expansion of the product about is

I'm not sure what the imaginary part is about...if you ignore it, here's the comparison you get between computing the product vs using the real part of above:

x = logspace(0,5,100)';
y = besseli(1,x).*besselk(1,x);
z = 1./2./x - 3./16./x./x - 45./256./x./x./x;

loglog(x,y,'-')
hold on;
loglog(x,z,'--')

Anyway, maybe this strategy is not at all useful for your case, I don't know.

[David Goodmanson]

Hello sir,

I don't have the impression you have written 'too much' and in fact I have a couple more questiions.

Could you mention what the physical situation is?

Something must distinguish the regions. It seems that each region has its own value of p determined in some fashion, is that correct?

For boundary conditions, is it as simple as for example

C4n I() + c5n K() = C6n I() + c7n K()

C4n dI()/dy + c5n dK()/dy = const* [ C6n dI()/dy + c7n dK()/dy ]
or
C4n d(yI())/dy + c5n d(yK())/dy = const* [ C6n d(yI())/dy + c7n d(yK())/dy ]

[hosein Javan]

Hi Mr.Lee. thanks for having time on replying. unfortunately in my case, the product of I1(x)*K1(x) is of no use to me. but I was thinking at something else, what do you think if we can change the unit of radian to a larger unit so that the input of bessel function reduces. if I replace 2*pi by "k", what happens to my equations?

[hosein Javan]

Hello again Mr. Goodmanson. thanks for your patience. yes the boundary condition is somehow like that. in fact it states that and must be continuous at each y1,...,y5. where "k" and "m" are constants for each "n" and has been defined before. and the regions are distinguished by y1,...,y5 boundaries. p is the frequency across only "x" and not "y", so for each region it is the same.

Mr. Goodmanson what about this idea: the input of bessel function dimension is radian which is "dimensionless", however I was wondering if we change that unit, what change should I do to the equations?