2
$\begingroup$

I have an equation:

w1[u_]=-(1/2) u^2 ea e0 Sin[2 a[z]] + K (a''[z])

with numerical values,

K = 6.2*10^-12; ea = 10; e0 = 8.85*10^-12; L = 1;

nd=NDSolve[{w1[1]==0,a[0]==0,a[L]==0},a[z],{z,0,L}];

gr1 = Plot[a[z] /. nd[[1]], {z, 0, L}, PlotRange -> Automatic]

The NDSolve should give the following output for different u, for u equal to zero a is 90 and by increasing u, a decreases till 0. Boundary conditions: @ z=0 and L a is 0 (very known problem),

enter image description here

However, the NDSolve is unstable and not giving the exact solution for the above program. Any suggestions or recommendations.

$\endgroup$
12
  • $\begingroup$ what is ele[1]? And where is a[u] used? $\endgroup$
    – Nasser
    Jan 10, 2023 at 13:06
  • $\begingroup$ @Nasser made a mistake editing the question. $\endgroup$
    – a019
    Jan 10, 2023 at 13:08
  • 2
    $\begingroup$ Well, I doubt that NDSolve is making an error, This tells me you have something wrong with your ode or with some of the numerical values you are using. Are you sure you are using same ode and parameter values used to generate the plot you are showing? $\endgroup$
    – Nasser
    Jan 10, 2023 at 13:23
  • 1
    $\begingroup$ @bmf reader.elsevier.com/reader/sd/pii/… $\endgroup$
    – a019
    Jan 10, 2023 at 13:29
  • 1
    $\begingroup$ I had quick look at the paper. Notice it says for fig 4, it is using the normalized z* which is z/h (do not know what h is. Planck's constant?) and your U is E*h. I think the problem is scaling/normalization issue as the ode look OK. !Mathematica graphics May be if you use the normalized values, the plot will show OK since h=6.62607015 × 10-34 which will make big difference in scaling. $\endgroup$
    – Nasser
    Jan 10, 2023 at 13:40

1 Answer 1

6
$\begingroup$

This problem can be solved with using the Euler wavelets collocation method described in our paper, we have

w1[u_] = -u^2 Sin[2 a[z]]/kk/2 + (a''[z]); ea = 10; e0 = 
 8.85*10^-12; kk = 6.2*10^-12/ea/e0; L = 1;
UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
       n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; zcol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y; A = Array[as, {nn}];
X2[y_] := A . Psi[y]; X1[y_] = A . int1[y] + c1; 
X0[y_] = A . int2[y] + c1 y + c2; eq = 
 Table[-(1/2) u^2 Sin[2 X0[y]]/kk + X2[y] == 0, {y, 
   xcol}]; bc = {X0[0] == Pi/2, X0[1] == Pi/2};


var = Join[A, {c1, c2}]; Do[
 soln[j] = 
   FindRoot[Join[eq, bc] /. u -> j, 
    Table[{var[[i]], 1/10}, {i, Length[var]}]];, {j, 0, 7}]
Plot[Evaluate[Table[2 90 /Pi X0[y] /. soln[j], {j, 0, 7}]], {y, 0, 1},
  Frame -> True, PlotLegends -> Table[j, {j, 0, 7}]]

    

Figure 1

$\endgroup$
11
  • 1
    $\begingroup$ wonderful stuff! $\endgroup$
    – bmf
    Jan 10, 2023 at 15:52
  • $\begingroup$ @AlexTrounev Can the same method be extended to similar solution but more complex problem such as for two coupled equations and four boundary conditions? The end solution is similar but starting from 2 and stabilizing at 90. I am trying NDSolve , shooting method but after certain range the NDSolve with shooting method can not find the solution. eq1: w1[u_]= L/K3(1/4 ((𝐾2βˆ’πΎ3 ) Sin[4 πœƒ[z]] πœ™'[z]^2+2 Sin[2 πœƒ[z]] (e0 ea 𝑒'[z]^2+(𝐾3βˆ’πΎ1 ) πœƒ'[z]^2+𝐾2 πœ™'[z]^2 ))+(𝐾1 Cos[πœƒ[z]]^2+𝐾3 Sin[πœƒ[z]] ) (dπœƒ'[z])/d𝑧)=0. $\endgroup$
    – a019
    Jan 26, 2023 at 15:17
  • $\begingroup$ @Alex Trounev eq2. w2= L/K3(Cos[πœƒ[z]](βˆ’2(𝐾2+(𝐾2βˆ’πΎ3 ) Cos[2 πœƒ[z]] Sin[πœƒ[z]] πœƒ'[z] πœ™'[z]+Cos[πœƒ[z]](𝐾2 Cos[πœƒ[z]]^2+𝐾3 Sin[πœƒ[z]]^2) (dπœ™'[z])/d𝑧))=0. The parameters are now, K1=13.4*10^-12, K2=6*10^-12, K3=18.5*10^-12, ea=11.35, e0=8.85*10^-12, L=1. Boundary conditions, πœƒ[0]== πœƒ[L]==2*(Pi/180), πœ™[0]==0, πœ™[L]==Pi/2. And u goes again from 0,to 7. The end solution is similar but starting from 2 and stabilizing at 90. Kind of mirror image of the above. Instead of a[z], now there is πœƒ[z] and πœ™[z] and {z,0,L}. Extension of the same problem. $\endgroup$
    – a019
    Jan 26, 2023 at 15:20
  • $\begingroup$ @MuhammadAli Yes, it can be extended to solve more complicated problem. $\endgroup$ Jan 27, 2023 at 4:26
  • $\begingroup$ @AlexTrounev I am not getting the same output, the one you showed above. Is there anything to do with the version of Mathematica I am using (11.0). I am adding a link to the imgur.com/a/6ijQjDU I tried with different things and plotrange but still not. $\endgroup$
    – a019
    Jan 27, 2023 at 8:41

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.