How to calculate the integration that the integrand should be determined by the values of variables?

Guu is the integration I want to calculate. However, the integrand Tuu must be determined by the values of x and y because the code

kzz=kz0(real(kz0)>=0&imag(kz0)>=0)

As a result, there appear errors. How to solve this problem and obtain the value of Guu? Many thanks!

The codes are as follows.

    clear all;
    syms kz x y d
    m = 2;
    dd=2.106*(m+1);
    vh = 4;
    mu = 11;
    delta = 8;
    HBAR = 1.05457266e-34;
    ME = 9.1093897e-31;
    ELEC = 1.60217733e-19;
    Kh = 2.106;
    vKh = [0,0,0;Kh,0,0;-Kh,0,0;0,Kh,0;0,-Kh,0];
    kc = sqrt(2*ME*ELEC/HBAR^2)*1e-10;
    ku = kc*sqrt(mu+delta);
    kd = kc*sqrt(mu-delta);
    a3 = [pi/Kh,pi/Kh,sqrt(2)*pi/Kh];
    kuu =@(x,y)[-ku*sin(x)*cos(y), -ku*sin(x)*sin(y), kz];
    n=0:m;
    for p=1:5;
        for q=1:5;
        tuu(p,q)=(sum((kuu(x,y) + vKh(p,:)).^2)-ku^2)*(p==q)+ kc^2*vh*sum(exp(i*n*sum((vKh(q,:)-vKh(p,:)).*a3)))/(m+1)*(p~=q);
        end
    end
    dtuu=det(tuu);
    kz0=vpasolve(dtuu,kz);
    kz0=double(kz0);
    kzz=kz0(real(kz0)>=0&imag(kz0)>=0);
    tuu1=subs(tuu,kz,kzz(1));
    tuu2=subs(tuu,kz,kzz(2));
    tuu3=subs(tuu,kz,kzz(3));
    tuu4=subs(tuu,kz,kzz(4));
    tuu5=subs(tuu,kz,kzz(5));
    tuu11=double(tuu1);
    tuu22=double(tuu2);
    tuu33=double(tuu3);
    tuu44=double(tuu4);
    tuu55=double(tuu5);
    nuu1=null(tuu11);
    nuu2=null(tuu22);
    nuu3=null(tuu33);
    nuu4=null(tuu44);
    nuu5=null(tuu55);
    piuu=[nuu1,nuu2,nuu3,nuu4,nuu5];
    pei=[1;0;0;0;0];
    A=piuu\pei;
    psiuu1=A(1)*nuu1(1)*exp(i*kzz(1)*d)+A(2)*nuu2(1)*exp(i*kzz(2)*d)+A(3)*nuu3(1)*exp(i*kzz(3)*d)+A(4)*nuu4(1)*exp(i*kzz(4)*d)+A(5)*nuu5(1)*exp(i*kzz(5)*d);
    psiuu2=A(1)*nuu1(2)*exp(i*kzz(1)*d)+A(2)*nuu2(2)*exp(i*kzz(2)*d)+A(3)*nuu3(2)*exp(i*kzz(3)*d)+A(4)*nuu4(2)*exp(i*kzz(4)*d)+A(5)*nuu5(2)*exp(i*kzz(5)*d);
    psiuu3=A(1)*nuu1(3)*exp(i*kzz(1)*d)+A(2)*nuu2(3)*exp(i*kzz(2)*d)+A(3)*nuu3(3)*exp(i*kzz(3)*d)+A(4)*nuu4(3)*exp(i*kzz(4)*d)+A(5)*nuu5(3)*exp(i*kzz(5)*d);
    psiuu4=A(1)*nuu1(4)*exp(i*kzz(1)*d)+A(2)*nuu2(4)*exp(i*kzz(2)*d)+A(3)*nuu3(4)*exp(i*kzz(3)*d)+A(4)*nuu4(4)*exp(i*kzz(4)*d)+A(5)*nuu5(4)*exp(i*kzz(5)*d);
    psiuu5=A(1)*nuu1(5)*exp(i*kzz(1)*d)+A(2)*nuu2(5)*exp(i*kzz(2)*d)+A(3)*nuu3(5)*exp(i*kzz(3)*d)+A(4)*nuu4(5)*exp(i*kzz(4)*d)+A(5)*nuu5(5)*exp(i*kzz(5)*d);
    Tuux=imag(conj(psiuu1)*diff(psiuu1,d)+conj(psiuu2)*diff(psiuu2,d)+conj(psiuu3)*diff(psiuu3,d)+conj(psiuu4)*diff(psiuu4,d)+conj(psiuu5)*diff(psiuu5,d));
    Tuu=subs(Tuux,d,dd);
    Guu=integral2(ku*sin(x)*Tuu,0.01,pi/2,0,pi/4)

Answer 1

[Walter Roberson]

Your tuu involves kz, x, and y. det() of that will involve kz, x, and y. You ask to vpasolve() with respect to kz, so the result is going to be in terms of x and y. But you then expect to be able to double() that result that is in terms of the unresolved variables x and y.

You could ask to vpasolve() without specifying a variable in order to get a set of numeric solutions. The equations involve sin() and cos() so there would typically be multiple solutions -- possibly 10 of them for each x and y value. Any one particular numeric solution would be of questionable worth.

[Henan Fang]

@Walter Roberson Sorry, I do not know how to revise the codes, only change vapsolve(dtuu,kz) to vapsolve(dtuu)? Can you tell me in details? I also can not understand that how to realize the codes "kzz=kz0(real(kz0)>=0&imag(kz0)>=0)", i.e. choose the solutions satisfy the condition its real part and imag part both >=0.

[Walter Roberson]

kzz=kz0(real(kz0)>=0&imag(kz0)>=0)

would be fine to select real and imaginary parts both greater than or equal to zero, provided that kz0 is numeric or is symbolic with a constant value (no free variables.)

vpasolve(dtuu)

would get you one numeric solution

>> sol = vpasolve(dtuu)
sol = 
  struct with fields:

      kz: [1×1 sym]
       x: [1×1 sym]
       y: [1×1 sym]

   sol.x, sol.y, sol.kz

However, there is no certainty that the one numeric solution for kz will have real and imaginary components that are non-negative, and probably will not.

For your purposes, you appear to be wanting all possible solutions. That is a problem as there are an infinite number of solutions, 10 for each combination of x and y values:

>> KZ0 = solve(dtuu,kz)
KZ0 =
  root(z^5 + z^4*((7902054562218934111819248704405*cos(y)^2*sin(x)^2)/316912650057057350374175801344 + (7902054562218934111819248704405*sin(x)^2*sin(y)^2)/316912650057057350374175801344 - 63275079194710162523/8796093022208000000) + z^3*(- (161656798990545805060539286473384785831975944601291*cos(y)^2*sin(x)^2)/696898287454081973172991196020261297061888000000 + (12488493260857014087641071366686064557292957871515854213280805*cos(y)^4*sin(x)^4)/50216813883093446110686315385661331328818843555712276103168 - (161656798990545805060539286473384785831975944601291*sin(x)^2*sin(y)^2)/696898287454081973172991196020261297061888000000 + (12488493260857014087641071366686064557292957871515854213280805*sin(x)^4*sin(y)^4)/50216813883093446110686315385661331328818843555712276103168 + (12488493260857014087641071366686064557292957871515854213280805*cos(y)^2*sin(x)^4*sin(y)^2)/25108406941546723055343157692830665664409421777856138051584 + 9202961831497445285815302223770634025269880912066396040017549882025405551/1606938044258990275541962092341162602522202993782792835301376000000000000) + z^2*((63606018950437609483403762545755898170282687640357634114174273629795398639697535926028871726917393305977*cos(y)^2*sin(x)^2)/101851798816724304313422284420468908052573419683296812531807022467719064988166835309169868800000000000 - (1058778775110342957793085448484573960427626989750064554192072757529036739296701617*cos(y)^4*sin(x)^4)/441711766194596082395824375185729628956870974218904739530401550323154944000000 + (19736951029439116276193759842601809558958615503769781801294745506428451490737303008141089205*cos(y)^6*sin(x)^6)/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792 + (63606018950437609483403762545755898170282687640357634114174273629795398639697535926028871726917393305977*sin(x)^2*sin(y)^2)/101851798816724304313422284420468908052573419683296812531807022467719064988166835309169868800000000000 - (1058778775110342957793085448484573960427626989750064554192072757529036739296701617*sin(x)^4*sin(y)^4)/441711766194596082395824375185729628956870974218904739530401550323154944000000 + (19736951029439116276193759842601809558958615503769781801294745506428451490737303008141089205*sin(x)^6*sin(y)^6)/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792 - (1058778775110342957793085448484573960427626989750064554192072757529036739296701617*cos(y)^2*sin(x)^4*sin(y)^2)/220855883097298041197912187592864814478435487109452369765200775161577472000000 + (59210853088317348828581279527805428676875846511309345403884236519285354472211909024423267615*cos(y)^2*sin(x)^6*sin(y)^4)/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792 + (59210853088317348828581279527805428676875846511309345403884236519285354472211909024423267615*cos(y)^4*sin(x)^6*sin(y)^2)/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792 + 4761821937671347871113761978335359455575618115183860939561896653203337937653498769553254268311/113078212145816597093331040047546785012958969400039613319782796882727665664000000000000000000) + z*((456540040115393944393393948615200500681044317192903792966196521548366841036733102883306648815649008981901228764612390378183*cos(y)^2*sin(x)^2)/17917957937422433684459538244547554224973163977877196279199912807710334969441287563047019946172856926208000000000000000000 + (936277788376230880901911312650896909654302647697483094487570578517551710857548450395777295603231171633969811784274792758373505764208237*cos(y)^4*sin(x)^4)/161390617380431786853494948250188242145606612051826469551916209783790476376052574664352834580008614464743948248296718336000000000000 - (711767051079138406303835086136470606414425388203685741613746534644374512112233507037930098518035557139732765867*cos(y)^6*sin(x)^6)/69992023193056381579920071267763883691301421788582797965624659405118495974380029543152421664737722368000000 + (31192492785294211383039122254098213894562309721869417651451953155101665863899566239792903548652224937058373650728356289605*cos(y)^8*sin(x)^8)/10086913586276986678343434265636765134100413253239154346994763111486904773503285916522052161250538404046496765518544896 + (456540040115393944393393948615200500681044317192903792966196521548366841036733102883306648815649008981901228764612390378183*sin(x)^2*sin(y)^2)/17917957937422433684459538244547554224973163977877196279199912807710334969441287563047019946172856926208000000000000000000 + (936277788376230880901911312650896909654302647697483094487570578517551710857548450395777295603231171633969811784274792758373505764208237*sin(x)^4*sin(y)^4)/161390617380431786853494948250188242145606612051826469551916209783790476376052574664352834580008614464743948248296718336000000000000 - (711767051079138406303835086136470606414425388203685741613746534644374512112233507037930098518035557139732765867*sin(x)^6*sin(y)^6)/69992023193056381579920071267763883691301421788582797965624659405118495974380029543152421664737722368000000 + (31192492785294211383039122254098213894562309721869417651451953155101665863899566239792903548652224937058373650728356289605*sin(x)^8*sin(y)^8)/10086913586276986678343434265636765134100413253239154346994763111486904773503285916522052161250538404046496765518544896 + (313581522049989043114180886947345407492978875308568008215413585373274964521989150796238082445132459831164679496028042924554289369745737*cos(y)^2*sin(x)^4*sin(y)^2)/16139061738043178685349494825018824214560661205182646955191620978379047637605257466435283458000861446474394824829671833600000000000 - (2135301153237415218911505258409411819243276164611057224841239603933123536336700521113790295554106671419198297601*cos(y)^2*sin(x)^6*sin(y)^4)/69992023193056381579920071267763883691301421788582797965624659405118495974380029543152421664737722368000000 - (2135301153237415218911505258409411819243276164611057224841239603933123536336700521113790295554106671419198297601*cos(y)^4*sin(x)^6*sin(y)^2)/69992023193056381579920071267763883691301421788582797965624659405118495974380029543152421664737722368000000 + (31192492785294211383039122254098213894562309721869417651451953155101665863899566239792903548652224937058373650728356289605*cos(y)^2*sin(x)^8*sin(y)^6)/2521728396569246669585858566409191283525103313309788586748690777871726193375821479130513040312634601011624191379636224 + (93577478355882634149117366762294641683686929165608252954355859465304997591698698719378710645956674811175120952185068868815*cos(y)^4*sin(x)^8*sin(y)^4)/5043456793138493339171717132818382567050206626619577173497381555743452386751642958261026080625269202023248382759272448 + (31192492785294211383039122254098213894562309721869417651451953155101665863899566239792903548652224937058373650728356289605*cos(y)^6*sin(x)^8*sin(y)^2)/2521728396569246669585858566409191283525103313309788586748690777871726193375821479130513040312634601011624191379636224 + (47667373761198201072851205569990915091189013139750604334500063275481*cos(y)*sin(x)^2*sin(y))/4074071952668972172536891376818756322102936787331872501272280898708762599526673412366794752000000 - 745168146905454103774898472131560592665009620057248732754445802686452517243692081012678857922334062215291437862139/7957171782556586274486115970349133441607298412757563479047423630290551952200534008528896000000000000000000000000) + (75334037659734271502216250392896192782195805548273506756897755552540204116567626982843865970638761*cos(y)^3*sin(x)^4*sin(y))/1291124939043454294827959586001505937164852896414611756415329678270323811008420597314822676640068915717951585986373746688000000 + (75334037659734271502216250392896192782195805548273506756897755552540204116567626982843865970638761*cos(y)*sin(x)^4*sin(y)^3)/1291124939043454294827959586001505937164852896414611756415329678270323811008420597314822676640068915717951585986373746688000000 - (2970920651603479611122317159822781144579257181170836469881522773157339216977793973761635371080759839373517471374985999372601475849123397049017254168781*cos(y)^2*sin(x)^4*sin(y)^2)/72683872429560689054932380788800453435364136068731806028149019918063928811339792332619105071376356556076252160626617793353460162861465600000000000 + (832306650974572558004690116779446632713548434225214129040734254956998507909604390126851849645643184244080393701989612139682314307082471644186595649731731984931574297*cos(y)^4*sin(x)^6*sin(y)^2)/10229345649675443343791217802586247350677006393884577467135285525300418113764607984010219038518450491096520887898625221903803926705891853291651648716800000000000 + (832306650974572558004690116779446632713548434225214129040734254956998507909604390126851849645643184244080393701989612139682314307082471644186595649731731984931574297*cos(y)^2*sin(x)^6*sin(y)^4)/10229345649675443343791217802586247350677006393884577467135285525300418113764607984010219038518450491096520887898625221903803926705891853291651648716800000000000 - (17529654853904749331101889792176989261990395063878524562885562532457086162637954953909*cos(y)*sin(x)^2*sin(y))/293567822846729153486185074598667128421960318613539983838411371441526128139326055432962374798096087878991872000000 - (262706954772456632468796746326776441473689206291692249888311144592769418060603404569344451552970625766682111727997482963458079700700531465*cos(y)^4*sin(x)^8*sin(y)^4)/2839213766779714416208296124562517712318911565184836172974571090549372219192960637992933791850638927971728600024477257552869537611776 - (87568984924152210822932248775592147157896402097230749962770381530923139353534468189781483850990208588894037242665827654486026566900177155*cos(y)^6*sin(x)^8*sin(y)^2)/1419606883389857208104148062281258856159455782592418086487285545274686109596480318996466895925319463985864300012238628776434768805888 - (87568984924152210822932248775592147157896402097230749962770381530923139353534468189781483850990208588894037242665827654486026566900177155*cos(y)^2*sin(x)^8*sin(y)^6)/1419606883389857208104148062281258856159455782592418086487285545274686109596480318996466895925319463985864300012238628776434768805888 + (49296955984203062054982667355315771504560039838949478973749622415775121968479383466986477944980981976171403653910394118993665179996305018201825843842005*cos(y)^8*sin(x)^10*sin(y)^2)/3196670515523576044934755563308202297086564498088930458479776726656380660551439995003193449537015778467662777468320381844938727095591204153641140224 + (49296955984203062054982667355315771504560039838949478973749622415775121968479383466986477944980981976171403653910394118993665179996305018201825843842005*cos(y)^6*sin(x)^10*sin(y)^4)/1598335257761788022467377781654101148543282249044465229239888363328190330275719997501596724768507889233831388734160190922469363547795602076820570112 + (49296955984203062054982667355315771504560039838949478973749622415775121968479383466986477944980981976171403653910394118993665179996305018201825843842005*cos(y)^4*sin(x)^10*sin(y)^6)/1598335257761788022467377781654101148543282249044465229239888363328190330275719997501596724768507889233831388734160190922469363547795602076820570112 + (49296955984203062054982667355315771504560039838949478973749622415775121968479383466986477944980981976171403653910394118993665179996305018201825843842005*cos(y)^2*sin(x)^10*sin(y)^8)/3196670515523576044934755563308202297086564498088930458479776726656380660551439995003193449537015778467662777468320381844938727095591204153641140224 - (87568984924152210822932248775592147157896402097230749962770381530923139353534468189781483850990208588894037242665827654486026566900177155*sin(x)^8*sin(y)^8)/5678427533559428832416592249125035424637823130369672345949142181098744438385921275985867583701277855943457200048954515105739075223552 - (87568984924152210822932248775592147157896402097230749962770381530923139353534468189781483850990208588894037242665827654486026566900177155*cos(y)^8*sin(x)^8)/5678427533559428832416592249125035424637823130369672345949142181098744438385921275985867583701277855943457200048954515105739075223552 - (668835488823690094606272667089432908390335339926913186919325473504982115756244248465086686679625596583880102887374170894633995578704346280634422881281*sin(x)^4*sin(y)^4)/726838724295606890549323807888004534353641360687318060281490199180639288113397923326191050713763565560762521606266177933534601628614656000000000000 - (668835488823690094606272667089432908390335339926913186919325473504982115756244248465086686679625596583880102887374170894633995578704346280634422881281*cos(y)^4*sin(x)^4)/726838724295606890549323807888004534353641360687318060281490199180639288113397923326191050713763565560762521606266177933534601628614656000000000000 + (9859391196840612410996533471063154300912007967789895794749924483155024393695876693397295588996196395234280730782078823798733035999261003640365168768401*sin(x)^10*sin(y)^10)/3196670515523576044934755563308202297086564498088930458479776726656380660551439995003193449537015778467662777468320381844938727095591204153641140224 + (9859391196840612410996533471063154300912007967789895794749924483155024393695876693397295588996196395234280730782078823798733035999261003640365168768401*cos(y)^10*sin(x)^10)/3196670515523576044934755563308202297086564498088930458479776726656380660551439995003193449537015778467662777468320381844938727095591204153641140224 - (2152931715275052911540607342144307993710121650723909911675232344873419884504941614723389769746799099574380152308566148766815064364152592678404459*sin(x)^2*sin(y)^2)/2521728396569246669585858566409191283525103313309788586748690777871726193375821479130513040312634601011624191379636224000000000000000000000000 - (2152931715275052911540607342144307993710121650723909911675232344873419884504941614723389769746799099574380152308566148766815064364152592678404459*cos(y)^2*sin(x)^2)/2521728396569246669585858566409191283525103313309788586748690777871726193375821479130513040312634601011624191379636224000000000000000000000000 + (721687976214858655313632357874108280699333481393247488098239494864472890976689140995908585282757110255400213482460577563240017625496739845657604886480269891588600703*sin(x)^6*sin(y)^6)/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584000000000000 + (721687976214858655313632357874108280699333481393247488098239494864472890976689140995908585282757110255400213482460577563240017625496739845657604886480269891588600703*cos(y)^6*sin(x)^6)/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584000000000000 + 1979287024797433578492127225054347055568880898653836667638803888767935838271523522314035102660950410833950100789924158366604738957/35835915874844867368919076489095108449946327955754392558399825615420669938882575126094039892345713852416000000000000000000000000, z, 1)^(1/2)

along with the negatives of those, and varying the root selector from 1 to 5, for a total of 10 different roots for each x and y.

Solutions for various values of x and y:

>> vpa(solve(simplify(subs(dtuu,{x,y}, {rand, rand}))))
ans =
   1.3140550680619409510293779556917
    1.604345239916760920828600880507
   1.9880575669901069090071341504563
  2.3925127528599078951466495854068i
  2.9212014474739945311628595705859i
  -1.3140550680619409510293779556917
   -1.604345239916760920828600880507
  -1.9880575669901069090071341504563
 -2.3925127528599078951466495854068i
 -2.9212014474739945311628595705859i
>> vpa(solve(simplify(subs(dtuu,{x,y}, {rand, rand}))))
ans =
  0.82736365817310165166497552876077
   1.2918420963741880296439925813063
    1.539152770917742652338779034751
  1.7033734463088766082621822940523i
   2.2299593435612285512014606421801
 -0.82736365817310165166497552876077
  -1.2918420963741880296439925813063
   -1.539152770917742652338779034751
 -1.7033734463088766082621822940523i
  -2.2299593435612285512014606421801

Solutions are periodic:

>> vpa(solve(simplify(subs(dtuu,{x,y}, {1/2, 3/4}))))
ans =
   1.4501422028519285986415170683356
  1.6514894522568748901797417494265i
   1.9099050759432311832908069854579
   1.9861673280320519970959561986972
  2.3464255011154026523172302884611i
  -1.4501422028519285986415170683356
 -1.6514894522568748901797417494265i
  -1.9099050759432311832908069854579
  -1.9861673280320519970959561986972
 -2.3464255011154026523172302884611i
>> vpa(solve(simplify(subs(dtuu,{x,y}, {1/2+2*pi, 3/4}))))
ans =
   1.4501422028519285654010393524238
  1.6514894522568741478008539691319i
   1.9099050759432310799461037354224
    1.986167328032052239422976502316
  2.3464255011154021802119485442698i
  -1.4501422028519285654010393524238
 -1.6514894522568741478008539691319i
  -1.9099050759432310799461037354224
   -1.986167328032052239422976502316
 -2.3464255011154021802119485442698i

You cannot list all of the possible values, only all of the possible values for any given x and y.

[Henan Fang]

@Walter Roberson I have changed vpasolve to solve, but there still appear errors. Can you teach me how to get Guu for my purpose in detail? Thank you very much.

[Walter Roberson]

You will have difficulty solving this using MATLAB. You can use solve() and index the first five values, which would often satisfy your range requirements, but I think that you would find it tricky to prove that the roots of the fifth order equation will definitely satisfy the requirements. The equation whose roots are being taken is long and complicated. Most of the trig terms are to even powers so most of the cos terms can be replaced by sqrt(1-sin^2) but there are a few with odd powers for which you would have to adjust the sign when you did that. Perhaps if you broke it up into the four sign cases for each of the two variables, a total of 16 different cases, then just maybe you might be able to prove that the roots must be in the right range... But I would tend to doubt it.

Just maybe you could establish that taking abs() of the roots would work. Maybe.

Anyhow, once roots were selected, they are going to be symbolic because of the x and y components. And as noted in one of your previous questions, null() does not work for symbolic expressions.

I think you probably have a definition problem in your equations. Because of the trig terms it is likely that for some x and y values that the null space is higher dimension than most of the cases. That would give you a piuu that was too large to \ against that particular vector.

[Henan Fang]

@Walter Roberson Thank you very much! I have voted your anwser.