witam pisze program w Maple, który metodą strzałów ma liczyć równanie schrodingera, cząstka umieszczona w kwadratowej studni potencjałów, niestety program w kilku miejscach i nie wysypuje, czy ktoś by pomógł?
> restart;with(plots): hbar:=6.6260755E-34; Me:=9.1093987E-31; ec:=1.60217733E-19; ck:=evalf((hbar*1E9/(2*Pi))^2 /(2*Me * ec));
-34
6.6260755 10
-31
9.1093987 10
-19
1.60217733 10
0.03809980339
>
> depth := 4; width := 2;
4
2
> V := -depth*(1-Heaviside(abs(x)-width/2)); g1:=plot(V,x=-width/2-1..width/2+1,color=blue,axes=framed,title=`V vs x`,labels=[nm,eV]): #`see if V looks ok` g1;
-4 + 4 Heaviside(|x| - 1)
> eq1 := -ck*(diff(y(x), x, x))+V*y(x) = eV*y(x);
/ d / d \\
-0.03809980339 |--- |--- y(x)||
\ dx \ dx //
+ (-4 + 4 Heaviside(|x| - 1)) y(x) = eV y(x)
> s_odd:= k -> dsolve({subs(eV=k,eq1),D(y)(0)=1,y(0)=0},y(x),numeric); s_even:= k -> dsolve({subs(eV=k,eq1),D(y)(0)=0,y(0)=1},y(x),numeric);
k -> dsolve({subs(eV = k, eq1), y(0) = 0, D(y)(0) = 1}, y(x),
numeric)
k -> dsolve({subs(eV = k, eq1), y(0) = 1, D(y)(0) = 0}, y(x),
numeric)
> test_point:=width/2+4; r_odd:= p -> (rhs(s_odd(p)(test_point)[2])); r_even:= p -> (rhs(s_even(p)(test_point)[2]));
5
p -> rhs(s_odd(p)(test_point)[2])
p -> rhs(s_even(p)(test_point)[2])
> _oddl:= p -> 3*sign(rhs(s_odd(p)(test_point)[2])); r_evenl:= p -> 3*sign(rhs(s_even(p)(test_point)[2]));
p -> 3 sign(rhs(s_odd(p)(test_point)[2]))
p -> 3 sign(rhs(s_even(p)(test_point)[2]))
> plot(r_oddl,-depth..0,y=-2..2,color=blue,labels=[eV,` `],ytickmarks=0,title=`odd energy levels`); plot(r_evenl,-depth..0,y=-2..2,color=blue,labels=[eV,` `],ytickmarks=0,title=`even energy levels`);
Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
> levels := [even, -3.9], [odd, -3.7], [even, -3.3], [odd, -2.8], [even, -2.1], [odd, -1.3], [even, -.4];
[even, -3.9], [odd, -3.7], [even, -3.3], [odd, -2.8],
[even, -2.1], [odd, -1.3], [even, -0.4]
> numlevels := nops([levels]);
7
> dr := proc (x) options operator, arrow; (r(x+0.1e-5)-r(x))/0.1e-5 end proc;
r(x + 0.000001) - r(x)
x -> ----------------------
0.000001
> for j to numlevels do end do;
> level:=levels[j][1]; mid:=levels[j][2];
Error, invalid subscript selector
Error, invalid subscript selector
> if level='odd' then r:=r_odd; s:=s_odd; else r:=r_even; s:=s_even; fi;
r_even
s_even
> mid0:=0: counter:=0: while (abs(mid-mid0)>0.000000002) and (counter < 10) do mid0:=mid: mid:=mid0-r(mid0)/dr(mid0): print(mid); counter:=counter + 1: od:
Error, cannot determine if this expression is true or false: 0.2e-8 < abs(mid)
> energy:=convert(mid,string); g2:=odeplot(s(mid),[x,y(x)],-width/2-1..width/2+1,axes=framed,title=`Ev = `.energy,labels=[nm,eV]): print(display([g2,g1]));
"mid"
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up
Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up
> print(s(mid.j)); mid.j:=mid; s.j:=s(mid.j); print(Ev=mid.j); od:
Error, unable to parse
printApplyFunction(sApplyFunction(midperiodj))semi midperiodj
Assignmidsemi speriodjAssignsApplyFunction(midperiodj)semi
printApplyFunction(Evequalsmidperiodj)semi
Typesetting:-mambiguous(od, Typesetting:-merror(unable to parse)
)colon
> for j from 1 to numlevels do s_abs.j:=(x) -> rhs(s.j(x)[2])^2; energy:=convert(mid.j,string); plot(s_abs.j,-width/2-1/2..width/2+1/2,title=`original |psi|^2 Ev=`.energy); area.j:=evalf(Int(s_abs.j,x=-width/2-1/2..width/2+1/2)); od;
2
x -> rhs(Typesetting:-delayDotProduct(s, j(x)[2]))
"mid"
Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
Error, (in evalf/int) invalid arguments
> for j from 1 to numlevels do s_norm.j:= (x) -> rhs(s.j(x)[2])/sqrt(area.j); energy:=convert(mid.j,string): g2:=plot(s_norm.j,evalf(-width/2-1..width/2+1),labels=[nm,eV],title=`normalized psi Ev=`.energy,axe s=framed): print(display([g2,g1])); od:
Error, (in plot) unexpected option: axe*s_even = framed