# Help to plot a surface plot with infinite roots?

**URL:** <https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291>\
**Category:** Visualization\
**Tags:** plotting, visualization, roots\
**Created:** [May 4, 2023, 3:00am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291 "2023-05-04T03:00:38Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Moslem\_Uddin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moslem_uddin/32/49483_2.png) [@Moslem\_Uddin](https://discourse.julialang.org/u/Moslem_Uddin)\
**Post date:** [May 4, 2023, 3:00am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/1 "2023-05-04T03:00:38Z")

</div>

I am trying to generate a surface plot of the following in Julia. As I’m new, I’m struggling  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/8/08d8fdf9cc16e548a46df667f17f8854c68948ae.png)  
where \alpha\_k are

I’m being able to generate a line plot at a fixed time using the following code:

```julia
using Roots,Plots
plotlyjs()

B= 1
Lambda_s=0
f(x) = tan(x) + x ./ (B*(1 .- Lambda_s * x .^2))

N=20
j0=zeros(N)
    
j0[1]=find_zero(f, 0)

# # print(f.(z0))

for i = 1:length(j0)-1
    j0[i+1] = j0[i] + 2.5; 
    j0[i+1] = find_zero(f,j0[i+1]);
end

# print(j0)
# print(f.(j0))
y=f.(j0)
 
xx = range(0,10,length=100) 
yy=f.(xx)
p1=plot(xx,yy, ylims=[-6, 6], title="Eigenvalues" , linewidth=2);
p2=plot(y,xlims=[0, 10],title="Function evaluated at the zeros" , linewidth=2);
 
p3=plot(j0,title="location of zeros" , linewidth=2);
plot(p1,p2,p3, legend=false)

t=0.2
x=range(0,1,length(j0))

f1= 1 .+ (cos.(j0) .^2 ./(B*(1 .- Lambda_s* j0 .^2))) .+ 2*B*Lambda_s*sin.(j0) .^2
fun= 2*((1 ./ j0) .* (sin.(j0 .* (1 .- x)) ./ f1)) .*exp.(-j0.^2*t)
ff= (1 .+ B*x)/(B .+ 1) .- sum(fun,dims=2)

plot(x,ff, linewidth=2)

```

Now, I would like to generate a surface plot between t=0 to t=2. Your help will be appreciated.

The above equations are from [https://link.springer.com/article/10.1007/s11012-020-01185-3](https://link.springer.com/article/10.1007/s11012-020-01185-3)

---

<div class="post-metadata">

**Author:** ![Moslem\_Uddin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moslem_uddin/32/49483_2.png) [@Moslem\_Uddin](https://discourse.julialang.org/u/Moslem_Uddin)\
**Post date:** [May 4, 2023, 3:01am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/2 "2023-05-04T03:01:38Z")

</div>

\alpha\_k are given by

![Screenshot 2023-05-03 215624](https://global.discourse-cdn.com/julialang/original/3X/3/9/3969cf7c2ee688e7f3dd830e760fefe2a3f1ff71.png)

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [May 4, 2023, 4:58am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/3 "2023-05-04T04:58:05Z")

</div>

Write the function just as it’s written in math (simplifying because on phone)

```julia
v(y, t) = A*y + sum(x->sin(B*x^2)*exp(-a[x]*t), 1:N)

```

And then

```julia
plot(v, (0, 10), (0, 2); seriestype=:surface)

```

If the series doesn’t converge nicely (for a small enough `N`), you’ll need to solve this further, possibly using an analytic framework. But this is how you’d do the plot part.

---

<div class="post-metadata">

**Author:** ![Moslem\_Uddin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moslem_uddin/32/49483_2.png) [@Moslem\_Uddin](https://discourse.julialang.org/u/Moslem_Uddin)\
**Post date:** [May 6, 2023, 5:47am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/4 "2023-05-06T05:47:25Z")

</div>

I tried the following:

```julia
# t=0.2
xs=range(0,1,length(j0))
ts=range(0,2,length(j0))
# analytic_sol_func(t,x) = (x )-(2/pi)*sum([(1/j0[k]* sin(pi*k*(1-x) )*exp(-k^2*pi^2*t ) for k in 1:1:500000])

# f1= 1 .+ (cos.(j0) .^2 ./(B*(1 .- Lambda_s* j0 .^2))) .+ 2*B*Lambda_s*sin.(j0) .^2
# fun= 2*((1 ./ j0) .* (sin.(j0 .* (1 .- x)) ./ f1)) .*exp.(-j0.^2*t)
# ff= (1 .+ B*x)/(B .+ 1) .- sum(fun,dims=2)
# sum([j0[k] for k in 1:1:length(j0)])
# fun= 2*((1 ./ j0) .* (sin.(j0 .* (1 .- x)) ./ (1 .+ (cos.(j0) .^2 ./(B*(1 .- Lambda_s* j0 .^2))) .+ 2*B*Lambda_s*sin.(j0) .^2))) .*exp.(-j0.^2*t)
ff(t,x)= (1 .+ B*x)/(B .+ 1) .- sum([2*((1 ./ j0[k]) .* (sin.(j0[k] .* (1 .- x)) ./ (1 .+ (cos.(j0[k]) .^2 ./(B*(1 .- Lambda_s* j0[k] .^2))) .+ 2*B*Lambda_s*sin.(j0[k]) .^2))).*exp.(-j0[k].^2*t) for k in 1:1:length(j0)])
u_real = reshape([ff(t,x) for t in ts for x in xs], (length(ts),length(xs)))

# surface(ts, xs, ff)
plot(u_real ,ts, xs; seriestype=:surface)

```

However, it seems like not working. Could you please explain a little more detailed way?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 6, 2023, 9:42am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/5 "2023-05-06T09:42:41Z")

</div>

When finding the roots of: `tan(x) + x = 0`, moving forward by steps of `2.5` is not enough to jump across the singularities. Try:

```julia
j0[i+1] = j0[i] + π

```

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [May 6, 2023, 8:34pm UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/6 "2023-05-06T20:34:11Z")

</div>

For future reference, it’s really difficult to debug something based on “it’s not working” and separate pieces of code that I don’t know how to stitch together. [Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757).

One thing I see immediately is that the surface plot takes the arguments in the order `x, y, Z`. But when I run your code, `u_real` is all `NaN`, so it’s probably not going to be a very pretty plot.

Here’s a rewrite from the math:

```julia
using Roots, Plots
B = 1
Λ_s = 0
N = 20

v(y, t; B, α, Λ_s) = (B*y + 1)/(B + 1) - 2*sum(α_k -> sin(α_k*(1-y))/(α_k*(1 + cos(α_k)^2/(B*(1-Λ_s*α_k^2)) + 2*B*Λ_s*sin(α_k)^2))*exp(-α_k^2*t), α)
f(x) = tan(x) + x/(B*(1-Λ_s*x^2))
α = zeros(N)
for k in 1:N
   α[k] = find_zero(f, (nextfloat(π*(k-3//2)), prevfloat(π*(k-1//2))))
end
plot(α)
plot(range(0,2,100), range(0,1,100), (y,t)->v(y, t; B, α, Λ_s); seriestype=:surface)

```

---

<div class="post-metadata">

**Author:** ![Moslem\_Uddin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moslem_uddin/32/49483_2.png) [@Moslem\_Uddin](https://discourse.julialang.org/u/Moslem_Uddin)\
**Post date:** [May 6, 2023, 11:11pm UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/7 "2023-05-06T23:11:29Z")

</div>

Thank you very much for your suggestion. However, when I tried your code, it gives me the following error:

GKS: Rectangle definition is invalid in routine SET\_WINDOW  
GKS: Rectangle definition is invalid in routine CELLARRAY  
invalid range

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [May 7, 2023, 3:41pm UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/8 "2023-05-07T15:41:10Z")

</div>

That’s because the function always returns `NaN`, because the first term in the summation has `1/α_k = 1/0`. Of course analytically we can recognize this as `sin(α)/α = 1`, but Julia cannot. We need to enforce that the first term is 0, for instance

```julia
v(y, t; B, α, Λ_s) = (B*y + 1)/(B + 1) - 2*sum(α_k -> α_k ≈ 0 ? 1 : sin(α_k*(1-y))/(α_k*(1 + cos(α_k)^2/(B*(1-Λ_s*α_k^2)) + 2*B*Λ_s*sin(α_k)^2))*exp(-α_k^2*t), α)

```

![v](https://global.discourse-cdn.com/julialang/original/3X/2/0/205c7632f6c78a9d32596c3db4da657edfa2b8e0.png)

Is that approximately what you expect?

Note that your quoted expression is defined for _all_ `α`, not just the `20` first positive `α` starting at `0.0`. I don’t know how well this converges.

---

<div class="post-metadata">

**Author:** ![Moslem\_Uddin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moslem_uddin/32/49483_2.png) [@Moslem\_Uddin](https://discourse.julialang.org/u/Moslem_Uddin)\
**Post date:** [May 15, 2023, 3:57am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/9 "2023-05-15T03:57:59Z")

</div>

Thanks for your suggestion.  
Avoiding zero gives me a better plot(((k+1/2)\*\pi\<\alpha[k]\<(k+3/2)\*\pi)

```julia
using Roots, Plots
plotlyjs()
B = 1
Λ_s = 1
N = 100

v(y, t; B, α, Λ_s) = (B*y + 1)/(B + 1) - 2*sum(α_k -> sin(α_k*(1-y))/(α_k*(1 + cos(α_k)^2/(B*(1-Λ_s*α_k^2)) + 2*B*Λ_s*sin(α_k)^2))*exp(-α_k^2*t), α)
# v(y, t; B, α, Λ_s) = (B*y + 1)/(B + 1) - 2*sum(α_k -> α_k ≈ 0 ? 1 : sin(α_k*(1-y))/(α_k*(1 + cos(α_k)^2/(B*(1-Λ_s*α_k^2)) + 2*B*Λ_s*sin(α_k)^2))*exp(-α_k^2*t), α)
f(x) = tan(x) + x/(B*(1-Λ_s*x^2))
α = zeros(N)

for k in 1:N
   α[k] = find_zero(f, (nextfloat(π*(k+3//2)), prevfloat(π*(k+1//2))))

end
# print(α)
plot(range(0,1,100), range(0,6,100), (y,t)->v(y, t; B, α, Λ_s); seriestype=:surface)

```

The above gives  
 ![newplot (11)](https://global.discourse-cdn.com/julialang/original/3X/e/6/e6d99d2ce1c3a8f152f8c8006d0ad33c0c965c51.png)  
While the true plot as reported in the article is  
 ![Screenshot 2023-05-14 211944](https://global.discourse-cdn.com/julialang/original/3X/f/1/f1eeaceb76092ba256dd5713af089819dde626a2.png)

Probably some other root-finding method will be needed. I would like to find roots in (0\<\alpha[k]\<(k+1/2)\*\pi). It seems like the above method is not working.

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [May 15, 2023, 6:08am UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/10 "2023-05-15T06:08:50Z")

</div>

Your limits appear to be in the wrong order, it should be `(lowerbound, upperbound)`. Do you know in which region they find their roots in the paper? Because I don’t know if this converges to the right thing if you only take the start of the positive axis. One strategy that might work is to take every other root negative (for instance, assign `\alpha[k]` and `\alpha[k+1]` in the loop, stepping `2:2:2N`.

---

<div class="post-metadata">

**Author:** ![Moslem\_Uddin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moslem_uddin/32/49483_2.png) [@Moslem\_Uddin](https://discourse.julialang.org/u/Moslem_Uddin)\
**Post date:** [May 21, 2023, 6:18pm UTC](https://discourse.julialang.org/t/help-to-plot-a-surface-plot-with-infinite-roots/98291/11 "2023-05-21T18:18:02Z")

</div>

It’s reported in article that, for 0\<\frac{1}{\sqrt{\Lambda\_s}}\leq \frac{\pi}{2}; \alpha\_k satisfy \pi k\<\alpha\_k\<(\pi+0.5)k. Motivated by this, I’ve tried the following:

```julia
using Roots, Plots
plotlyjs()
B = 1
Λ_s = 1
N = 100

v(y, t; B, α, Λ_s) = (B*y + 1)/(B + 1) - 2*sum(α_k -> sin(α_k*(1-y))/(α_k*(1 + cos(α_k)^2/(B*(1-Λ_s*α_k^2)) + 2*B*Λ_s*sin(α_k)^2))*exp(-α_k^2*t), α)

f(x) = tan(x) + x/(B*(1-Λ_s*x^2))
α = zeros(N)

for k in 1:N
   α[k] = find_zero(f, (nextfloat(π*(k)), prevfloat(π*(k+1//2))))
end
# print(α)
plot(range(0,1,N), range(0,6,N), (y,t)->v(y, t; B, α, Λ_s); seriestype=:surface)

```

I’m getting the following plot  
 ![newplot (11)](https://global.discourse-cdn.com/julialang/original/3X/e/6/e6d99d2ce1c3a8f152f8c8006d0ad33c0c965c51.png)  
It seems pretty reasonable, but not the desired one.