-
Notifications
You must be signed in to change notification settings - Fork 2
/
solutionM3.m
55 lines (35 loc) · 1.53 KB
/
solutionM3.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
clc
h=solve('x^3-(exp((-(3*512/10^9+512/10^9))*(r-r*x)))=0','x')
% format long
root=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];%用于存储z_m的根。
i=1;
rou=0.05;
while (rou<0.75)
r=rou*10^9/512;
root(i)=(-(5859375*exp(-(4*r)/1953125)^(1/3)*exp(- (pi*2i)/3 - log(exp(-(4*r)/1953125))/3)*((3^(1/2)*1i)/2 - 1/2)*lambertw(0, -(4*r*exp((pi*2i)/3 + log(exp(-(4*r)/1953125))/3))/5859375))/(4*r));
rou=rou+0.05;
i=i+1;
end
root
root1=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];%用于存储z_m的根。
i=1;
rou=0.05;
while (rou<0.75)
r=rou*10^9/512;
r;
root1(i)=((5859375*exp(-(4*r)/1953125)^(1/3)*exp(- (pi*4i)/3 - log(exp(-(4*r)/1953125))/3)*((3^(1/2)*1i)/2 + 1/2)*lambertw(0, -(4*r*exp((pi*4i)/3 + log(exp(-(4*r)/1953125))/3))/5859375))/(4*r));
rou=rou+0.05;
i=i+1;
end
root1
i=1;
while (i<14)
a=exp(-(512*10^(-9)*(i*0.05*10^9/512-(i*0.05*10^9/512)*root(i))));
b=exp(-512*10^(-9)*(i*0.05*10^9/512-(i*0.05*10^9/512)*root(i)));
a1=exp(-(512*10^(-9)*(i*0.05*10^9/512-(i*0.05*10^9/512)*root1(i))));
b1=exp(-512*10^(-9)*(i*0.05*10^9/512-(i*0.05*10^9/512)*root1(i)));
mat1=sym([3,2,1;((root(i))^3*(a^0)-a^3*(root(i))^0),((root(i))^3*a-a^3*root(i)),(root(i))^3*(a^2)-a^3*(root(i))^2;((root1(i))^3*(a1^0)-a1^3*(root1(i))^0),((root1(i))^3*a1-a1^3*root1(i)),(root1(i))^3*(a1^2)-a1^3*(root1(i))^2]);
mat2=sym([(3-(0.05*i*10^9/512)*(512/10^9)/(1-0.05*i));0;0]);
hmat=eval(mat1\mat2) %即方程组的解。
i=i+1;
end