0

I have such a code:

θ[x_] := Boole[x >= 0]
{T, Um} := {10^(-3), 10};
{t1, t2, t3} := {T/2, T/3, T/6};

u1[t_] := Um/t3(tθ[t] - (t - t3) θ[t - t3] - (t - t2) θ[t - t2] + (t - t1) θ[t - t1])
U1[s_] := LaplaceTransform[u1[t], t, s]
A1[ω_] := Abs[U1[s]] /. s -> I*ω


NCriterion := 1/10
A1max := Limit[A1[ω], ω -> 0]


ωPetal = Solve[A1[ω] == 0 && ω ∈ Interval[{10000, 20000}], ω][[1]];
ωPetal
ωNCriterion = Solve[A1[ω] == A1max*NCriterion && ω ∈ Interval[{10000, 20000}], ω][[1]];
ωNCriterion

When processing ωPetal everything is okay but when processing ωNCriterion I get a warning:

Solve was unable to prove that the solution set found is complete

and some strange output: Scary output. Why I observe such a behaviour of system and how do I make it process ωNCriterion as good as ωPetal? Thanks in advance.

Update: I need exact root value like {ω → 5246π}, not approximation like 16480.3, which I get with NSolve[] and FindRoot[].

Cpp Nosavvier
  • 67
  • 7
  • 1
    Perhaps an interesting comment: NSolve seems to be working fine producing a result equal to 16480.3. This s roughly equal to 5246 Pi. –  May 13 '21 at 12:13
  • 1
    Extra comment: FindRoot also works FindRoot[A1[[Omega]] - A1max*NCriterion == 0, {[Omega], 16000}] –  May 13 '21 at 12:22
  • @DiSp0sablE_H3r0 Thank you, but how to explain why Solve[] gives up and how do I get exact root value as {ω → 5246π} or something like that? – Cpp Nosavvier May 13 '21 at 13:42
  • 1
    Not really an expert to try and answer this one. Sorry. Perhaps something to do with internal structure of Solve or something related? Be that as it may, you can point out that NSolve and FindRoot give the same answer, and rephrase the question. Perhaps someone else might know :-) –  May 13 '21 at 13:56
  • @DiSp0sablE_H3r0 Okay, thanks again :) – Cpp Nosavvier May 13 '21 at 14:02
  • 1
    TheRoot[..] is the exact answer. See the answers here for more about Root[] objects. – Michael E2 May 13 '21 at 14:20
  • Further to @MichaelE2's comment: Not all solutions to algebraic equations can be expressed as rational numbers (or, in your case, rational multiples of π.) For example, the equation $x^5 - x + 1 = 0$ doesn't have solutions that can be expressed in terms of roots, fractions, etc. I suspect that your equation is one of these. – Michael Seifert May 13 '21 at 14:26

2 Answers2

3
Clear["Global`*"]

θ[x_] := Boole[x >= 0]
{T, Um} = {10^-3, 10};
{t1, t2, t3} = {T/2, T/3, T/6};

Using Set rather than SetDelayed in the definition of u1 results in the calculation being done once. And Simplify will produce a Piecewise expression.

u1[t_] = Um/
    t3*(t*θ[t] - (t - t3) θ[t - t3] - (t - t2) θ[
       t - t2] + (t - t1) θ[t - t1]) // Simplify

enter image description here

Similarly with U1

U1[s_] = LaplaceTransform[u1[t], t, s]

(* (60000 E^(-s/2000) (-1 + E^(s/6000))^2 (1 + E^(s/6000)))/s^2 *)


A1[ω_] = Abs[U1[s]] /. s -> I*ω;


NCriterion = 1/10;
A1max = Limit[A1[ω], ω -> 0];

Instead of ω ∈ Interval[{10000, 20000}] use either Between[ω, {10000, 20000}] or 10000 <= ω <= 20000

ωPetal = 
 Solve[A1[ω] == 0 && 10000 <= ω <= 20000, ω][[1]]

(* {ω -> 6000 π} *)

Note that the RHS of the Rule is no longer a List

ωNCriterion = Solve[A1[ω] == A1max*NCriterion && 
    10000 <= ω <= 20000, ω][[1]]

enter image description here

The Root expression is an exact solution. Its approximate numeric value is obtained with N

ωNCriterion // N

(* {ω -> 16480.3} *)


Plot[{A1[ω], A1max*NCriterion}, {ω, 10000, 20000},
 PlotLegends -> Placed["Expressions", {0.7, 0.6}],
 Epilog -> {Red, AbsolutePointSize[4], Tooltip @@ 
    N[{Point[{ω, A1[ω]}], {ω, A1[ω]}} /. ωNCriterion]}]

enter image description here

Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
1

Here's a somewhat simpler expression for the solution:

ωNCriterion /. HoldPattern@Root[{Function[body_], _}] :>
  (6000 ArcCos[x] /. 
    First@Solve[
      body == 0 && -1 < x < 0 /. {#1 -> 6000 ArcCos[x]} // Simplify, 
      x])
% // N
(*
{ω -> 
  6000 ArcCos[
    Root[{
     -ArcCos[#1]^2 + 10 Sqrt[2] Sqrt[(-1 + #1)^2 (1 + #1)] &, 
     -0.92304340472946770158}]
    ]}

{ω -> 16480.3}
*)
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thank you. eventually I preferred to solve the problem using ωNCriterion /. x_?NumericQ :> π*Round[N[x/π]]. – Cpp Nosavvier May 23 '21 at 17:06
  • @CppNosavvier In that case, FindRoot would be faster than Solve: Pi*Round[ω/Pi /. FindRoot[A1[ω] == A1max*NCriterion, {ω, 15000, 10000, 20000}]]. If you want a more accurate solution, you can take more of the terms of the continued fraction: Pi*FromContinuedFraction[Take[ContinuedFraction[ω/Pi /. FindRoot[A1[ω] == A1max*NCriterion, {ω, 15000, 10000, 20000}]], 3]] – Michael E2 May 23 '21 at 17:53