求解带有变量的积分方程组

Question

已知kB，va，hbar，T和kl，通过T和kl一一对应的关系，列出3个方程组，求解L,A,B的值

clear;clc;
%Sample=[MgAgSb];
vl=[3360]; %longitudinal sound velocities; 纵波声速，单位:m/s,列矩阵中间用分号;间隔
vt=[1715]; %transverse sound velocities; 纵波声速，单位:m/s
n=[3]; % number of atoms per formular ,化学式中原子数，无单位
p=[6.31]; % 密度 单位：g/cm3
M=[253.9332]; %化学式中相对分子质量，单位：g/mol
T=[4.37634,5.84135,7.25596]; % 温度 单位：K
klattice=[1.33336,1.93658,2.36888];
h=6.62607015e-34; %h presents the Planck constant.%单位(J·s)
hbar=1.05457266e-34; % reduced Planck constant 单位(
k=1.3806505e-23;%玻尔兹曼常数k，e在matlab中代表10，单位(J·K-1)
NA=6.02214076e23; % 阿伏伽德罗常数，单位：1/mol
va=((1/3).*((1./vl.^3)+(2./vt.^3))).^(-1/3); %平均声速
DebyeT=(h/k).*(va).*(((3e6/(4*(pi))).*(n.*NA.*p./M)).^(1/3)); % Debye temperature 单位:K

```bash
[L,A,B]=solve((k./(2*((pi)^2).*va)).*((k./hbar).^3).*(T(1).^3).*integral(@(x)((x.^4).*exp(x)./(((va./L)+(A.*(k.*T(1)./hbar).^4.*(x.^4))+B.*T(1).*exp(-DebyeT./(3*T(1))).*(k.*T(1)./hbar).^2.*x.^2).*((exp(x)-1).^2))),x,0,DebyeT/T(1))==klattice(1)',(k./(2*((pi)^2).*va)).*((k./hbar).^3).*(T(2).^3).*integral(@(x)((x.^4).*exp(x)./(((va./L)+(A.*(k.*T(2)./hbar).^4.*(x.^4))+B.*T(2).*exp(-DebyeT./(3*T(2))).*(k.*T(2)./hbar).^2.*x.^2).*((exp(x)-1).^2))),x,0,DebyeT/T(2))==klattice(2),(k./(2*((pi)^2).*va)).*((k./hbar).^3).*(T(3).^3).*integral(@(x)((x.^4).*exp(x)./(((va./L)+(A.*(k.*T(3)./hbar).^4.*(x.^4))+B.*T(3).*exp(-DebyeT./(3*T(3))).*(k.*T(3)./hbar).^2.*x.^2).*((exp(x)-1).^2))),x,0,DebyeT/T(3))==klattice(3),L,A,B); % 根据方程组求解值

运行结果及报错内容 ：

错误使用 integral (line 85)
A 和 B 必须为浮点标量。

出错 thermaltot (line 30)
[L,A,B]=solve((k./(2*((pi)^2).*va)).*((k./hbar).^3).*(T(1).^3).*integral(@(x)((x.^4).*exp(x)./(((va./L)+(A.*(k.*T(1)./hbar).^4.*(x.^4))+B.*T(1).*exp(-DebyeT./(3*T(1))).*(k.*T(1)./hbar).^2.*x.^2).*((exp(x)-1).^2))),x,0,DebyeT/T(1))==klattice(1)',(k./(2*((pi)^2).*va)).*((k./hbar).^3).*(T(2).^3).*integral(@(x)((x.^4).*exp(x)./(((va./L)+(A.*(k.*T(2)./hbar).^4.*(x.^4))+B.*T(2).*exp(-DebyeT./(3*T(2))).*(k.*T(2)./hbar).^2.*x.^2).*((exp(x)-1).^2))),x,0,DebyeT/T(2))==klattice(2),(k./(2*((pi)^2).*va)).*((k./hbar).^3).*(T(3).^3).*integral(@(x)((x.^4).*exp(x)./(((va./L)+(A.*(k.*T(3)./hbar).^4.*(x.^4))+B.*T(3).*exp(-DebyeT./(3*T(3))).*(k.*T(3)./hbar).^2.*x.^2).*((exp(x)-1).^2))),x,0,DebyeT/T(3))==klattice(3),L,A,B); 
>> 

之前尝试过quad积分和定义syms L,A,B,x为未知量带入也无法算算出来，


想要达到的结果：通过T和kl一一对应的关系，列出3个方程组，求解L,A,B的值
integral为啥无法积分出含有变量的定积分？

Answer 1

你好，既然有三个方程，就分开写，积分的时候上下限保持标量即可

Answer 2

积分上下限是标量，但是积分内部有L,A,B的变量

Answer 3

clear;clc;
vl=3360; %longitudinal sound velocities; 纵波声速，单位:m/s,列矩阵中间用分号;间隔
vt=1715; %transverse sound velocities; 纵波声速，单位:m/s
n=3; % number of atoms per formular ,化学式中原子数，无单位
p=6.31; % 密度 单位：g/cm3
M=253.9332; %化学式中相对分子质量，单位：g/mol
T1=4.376347; % 温度 单位：K
T2=5.84135;
T3=25596;
klattice1=1.33336;
klattice2=1.93658;
klattice3=2.36888;
h=6.62607015e-34; %h presents the Planck constant.%单位(J·s)
hbar=1.05457266e-34; % reduced Planck constant 单位(J·s)
k=1.3806505e-23;%玻尔兹曼常数k，e在matlab中代表10，单位(J·K-1)
NA=6.02214076e23; % 阿伏伽德罗常数，单位：1/mol
va=((1/3).((1./vl.^3)+(2./vt.^3))).^(-1/3); %平均声速
DebyeT=(h/k).(va).(((3e6/(4(pi))).(n.NA.p./M)).^(1/3)); % Debye temperature 单位:K
BM=p.(3(vl.^2)-4.(vt.^2))/3e6; %Bulk modulus 体积模量，单位：GPa
rG=(3/2).(3vl.^2-4vt.^2)./(vl.^2+2.vt.^2); % Gruneisen parameter
T11=DebyeT/T1;
T21=DebyeT/T2;
T31=DebyeT/T3;
syms x w e r
[w,e,r]=solve((k./(2((pi)^2).va)).((k./hbar).^3).(T1.^3).integral(@(x)((x.^4).exp(x)./(((va./w)+(e.(k.T1./hbar).^4.(x.^4))+r.T1.exp(-DebyeT./(3T1)).(k.T1./hbar).^2.x.^2).((exp(x)-1).^2))),x,0,T11)==klattice1,(k./(2((pi)^2).va)).((k./hbar).^3).(T2.^3).integral(@(x)((x.^4).exp(x)./(((va./w)+(e.(k.T2./hbar).^4.(x.^4))+r.T2.exp(-DebyeT./(3T2)).(k.T2./hbar).^2.x.^2).((exp(x)-1).^2))),x,0,T21)==klattice2,(k./(2((pi)^2).va)).((k./hbar).^3).(T3.^3).*integral(@(x)((x.^4).exp(x)./(((va./w)+(e.(k.T3./hbar).^4.(x.^4))+r.T3.exp(-DebyeT./(3T3)).(k.*T3./hbar).^2.x.^2).((exp(x)-1).^2))),x,0,T31)==klattice3,w,e,r);

我更新了代码，积分上下限都是标量，在上述代码中加粗表示了，但还是报错错误使用 integral (line 85) A 和 B 必须为浮点标量。
[w,e,r]=solve((k./(2*((pi)^2).va)).((k./hbar).^3).*(T1.^3).integral(@(x)((x.^4).exp(x)./(((va./w)+(e.(k.T1./hbar).^4.(x.^4))+r.T1.exp(-DebyeT./(3T1)).(k.T1./hbar).^2.x.^2).((exp(x)-1).^2))),x,0,46)==klattice1,(k./(2((pi)^2).va)).((k./hbar).^3).(T2.^3).integral(@(x)((x.^4).exp(x)./(((va./w)+(e.(k.T2./hbar).^4.(x.^4))+r.T2.exp(-DebyeT./(3T2)).(k.T2./hbar).^2.x.^2).((exp(x)-1).^2))),x,0,34)==klattice2,(k./(2((pi)^2).va)).((k./hbar).^3).(T3.^3).*integral(@(x)((x.^4).exp(x)./(((va./w)+(e.(k.T3./hbar).^4.(x.^4))+r.T3.exp(-DebyeT./(3T3)).(k.*T3./hbar).^2.x.^2).((exp(x)-1).^2))),x,0,0.0079)==klattice3,w,e,r);
甚至将积分上下限改为数值，也报错错误使用 integral (line 85)
A 和 B 必须为浮点标量。

请问是什么原因呀？

Answer 4

你把你的方程组写成类似这种

eq = @(x) [integral(f, x(1), x(2)) + x(3);
integral(f, x(2), x(3)) + x(1);
integral(f, x(3), x(1)) + x(2);];
x = fsolve(eq, x0)